NPS IMD Intro to R
NPS Inventory and Monitoring Division Intro to R Training: January 31 – February 2, 2022
Prep for Training
Installing required software
The only thing you really need to do to prepare for R training is to install the latest version of R and RStudio. We’ll talk about the difference between R and RStudio on the first day, but for now, just make sure they’re installed. Directions for installing R/RStudio are below. If you run into any problems, check the instructions at R for Data Science Section 1.4.1 and 1.4.2.
NOTE: If you don’t have Administrative privileges on your computer, you will have to submit an IT HelpDesk Ticket (need to be on VPN or NPS network to access this link) to install R and RStudio, which could take a while. PLEASE PLAN AHEAD!
Even if you already have R or RStudio installed, please install the latest versions of both programs. R recently went through a major version change from 3.x to 4.x with some potential code-breaking changes. The latest versions are needed to ensure everyone’s code behaves the same way.
Install the latest version of R (4.1.2 as of January 2022)
Assuming you’re on a PC, the installer for the latest version of R can be downloaded from The Comprehensive R Archive Network (CRAN). You’ll want to download the R installer by clicking on “Download R 4.0.4 for Windows”, or whatever the latest version number happens to be. After it downloads, open and run the installer with default settings. Note that there’s no need to pin R to your taskbar or add a shortcut to your desktop. For 99% of your needs, you’ll run R within RStudio, and won’t ever open R directly.Install RStudio
The installer for RStudio for PCs can be downloaded here. You’ll need to click on the large blue “DOWNLOAD RSTUDIO FOR WINDOWS” button to download the installer. After it downloads, open and run the installer with default settings. I like having RStudio pinned to my taskbar, so it’s easier to find/open, but it’s up to you and your preferences.Optional Reading
We know you’re all busy and have to prioritize how you spend your time. If you have any time before training starts, we highly recommend reading Chapter 2: R basics and workflows in STAT545. This is an online book based on a graduate level statistics class that Dr. Jenny Bryan teaches and is turning into this book. She’s a stats professor turned senior programmer at RStudio, and I really like how she approaches programming R in this book.Structure of training
The training will take place over 4 half days. Each day will run from 10-2 MST via MS Teams. The hour before training, and the hour after training will also be posted by at least one trainer as office hours, in case there are questions that couldn’t be handled during training.
Each training session has three trainers assigned, two of which will be the main instructors. The third trainer will manage the chat. For most of the training, one of the trainers will share their screen to walk through the website and then demo with live coding. It will help a lot if you have 2 monitors, so you can have the screen being shared on one monitor and your own session of RStudio open on the other.
Four days is barely enough to scratch the surface on what you can do in R. Our goals with this training are to help you get beyond the initial learning curve that can be really challenging to climb on your own, and to expose you to some of the really useful things R can do for you in your work. The trainers put lot of thought and time into designing this training. We did it because we enjoy coding in R and it has greatly increased efficiency and productivity in our jobs. We hope you have a similar experience as you start out on this R learning curve. As you learn more, please share your skills and code with others to help us build a larger community of R users in IMD who can learn from each other.
Finally, to help us improve this training for future sessions, please leave feedback in the training feedback form.
Day 1: Intro to R
Intro to R & RStudio
Some Notes on Coding
Welcome to the class!
Learning (or re-learning) how to program can feel like an intimidating learning curve. Congrats on taking the first step! Before we delve into technical topics, let’s talk a little about coding in general.
It’s okay to make mistakes
Errors and broken code are not only part of the learning process, they’re part of the coding process. If you think your code is running perfectly on the first try, you probably didn’t test it thoroughly enough. It can be frustrating and discouraging when you can’t get your code to work the way you want it to, but getting stuck doesn’t make you a bad programmer. Try to approach broken code as a puzzle and a learning opportunity.
Use the resources available to you
As the saying goes,
Good coders write good code; great coders steal great code.
Google and Stack Overflow are great resources when you don’t know off the top of your head how to do something. And reach out to your colleagues too - there’s no reason to waste time writing code that someone else has already written! Just give credit where credit is due, and take some time to make sure you understand what the copied code does. And especially don’t hesitate to reach out to more experienced colleagues with questions.
About R
R is a programming language that was originally developed by statisticians for statistical computing and graphics. R is free and open source, which means you will never need a paid license to use it, and anyone who wants to can view the underlying source code and suggest fixes and improvements. Since its first official release in 1995, R remains one of the leading programming languages for statistics and data visualization, and its capabilities continue to grow.
Every programming language has strengths and weaknesses. We like R because:
- It is great for statistics and data exploration
- It is good at producing publication-quality plots and tables
- There are many well-maintained add-ons (called packages) that add useful features to R (this website was built using the
rmarkdownpackage) - R is well-documented, with many learning resources available whether you are a beginner or an advanced user
For more information on the history of R, visit the R Project website.
About RStudio
RStudio is what’s called an integrated development environment (IDE). When you install R, it comes with a simple user interface that lets you write and execute code. Writing code in this interface is kind of like doing word processing in Notepad: it’s simple and straightforward, but soon you’ll start wishing for more features. This is where RStudio comes in. RStudio makes programming in R easier by color coding different types of code, autocompleting code, flagging mistakes (like spellcheck for code), and providing many useful tools with the push of a button or key stroke (e.g. viewing help info).
When you open RStudio, you typically see 4 panes:
Source:
This is where the code gets written. When you create a new script or open an existing one, it displays here. In the screenshot above, there’s a script called bat_data_wrangling.R open in the source pane. Note that if you haven’t yet opened or created a new script, you won’t see this pane until you do.
The source pane will helpfully color-code your code to make it easier to read. It will also detect syntax errors (the coding equivalent of a grammar checker). These are flagged with a red “x” next to the line number and a squiggly line under the offending code. When you’re ready to run all or part of your script:- Highlight the line(s) of code you want to run
- Either click the “Run” button (top right of the source pane) or press Ctrl+Enter. At this point, the code gets sent to the console (the bottom left pane). You’ll see your code appear in the console, and then you’ll see the output if there is any.
Console:
This is where the code actually runs. When you first open RStudio, the console will tell you the version of R that you’re running (should be R 4.1.0 or greater). We talked above about how you can run a script at the console. You can also type code directly into the console and run it that way. That code won’t get saved to a file, but it’s a great way to experiment and test out lines of code before you add them to your script in the source pane.
The console is also where errors appear if your code breaks. Deciphering errors can be a challenge sometimes, but Googling them is a good place to start.Environment/History/Connections:
- Environment: This is where you can see what is currently in your environment. Think of the environment as temporary storage for objects - things like datasets and stored values - that you are using in your script(s). You can also click on objects and view them.
Remember that anything you see in your environment is temporary and it will disappear when you restart R. If there is something in your environment that you want to access in the future, make sure that your script is able to reproduce it (or save it to a file). - History: This shows the code you’ve run in the current session. It’s not good to rely on it, but it can be a way to recover code you ran at the console and later realized you needed in your script.
- Connections: This is one way to connect R to a database. We’ll talk about this more later on.
- Git: If you have installed Git on your computer, you may see a Git tab. We won’t talk much about it this week, but this is where you’ll keep track of changes to your code. You are encouraged to attend the Git/GitHub session during Week 2 if you’d like to learn more!
- Tutorial: This has some interactive tutorials that you can check out if you are interested.
- Environment: This is where you can see what is currently in your environment. Think of the environment as temporary storage for objects - things like datasets and stored values - that you are using in your script(s). You can also click on objects and view them.
Files/Plots/Packages/Help/Viewer:
- Files: This tab shows the files within your working directory (typically the folder where your current code project lives). More on this later.
- Plots: This tab will show plots that you create.
- Packages: This tab allows you to install, load, and update packages, and also view the help files within each package. This is kind of like the R equivalent of an app store. You can also do this with code.
- Help: Allows you to search for and view documentation for packages that are installed on your computer.
- Viewer: Shows HTML outputs produced by R Markdown, R Shiny, and some plotting and mapping packages.
Global Options
There are several settings in the Global Options that everyone should check to make sure we all have consistent settings. Go to Tools -> Global Options and follow the steps below.
- Under the General tab, you should see that your R Version is [64-bit] and the version is R-4.1.0 or greater. If it’s not, you probably need to update R. Let an instructor know if you have questions about this.
- Make sure you are not saving your environment. To do this, uncheck the Restore .RData into your workspace at startup option.
When this option is checked, R Studio saves your current working environment (the stuff in the Environment tab) when you exit. The next time you open R Studio, it restores that same environment. This seems like a good thing, but the whole point of using R is that your code should return the same results every time you run it. Clearing your environment every time you close RStudio forces you to run your code with a clean slate so that you can be sure that it will work for everyone, not just you and your specific environment.
If you sometimes work with huge datasets that take a long time to load and process, you may want to set Save workspace to .RData on exit to “Ask”. When you close RStudio, it will ask you if you want to save your workspace image. If you have a huge dataset (or the products of a huge dataset) in your environment, select “Yes” and then runload(".RData")at the console the next time you open RStudio to restore everything.
Most other settings are whatever you prefer, including everything in the Appearance window.
Starting a Project
When you’re working on a code project it’s important to stay organized. Keeping code, input data, and results together in one place will protect your sanity and the sanity of the person who inherits the project. R Studio projects help with this. Creating a new R Studio project for each new code project makes it easier to manage settings and file paths.
Before we create a project, take a look at the Console tab. Notice that at the top of the console there is a folder path. That path is your current working directory.
If you refer to a file in R using a relative path, for example ./data/my_data_file.csv, R will look in your current working directory for a folder called data containing a file called my_data_file.csv. Using relative paths is a good idea because the full path is specific to your computer and likely won’t work on a different computer. But there’s no guarantee that everyone has the same default R working directory. This is where projects come in.
Let’s make a project for this class. Click File > New Project. In the window that appears, select New Directory, then select New Project. You will be prompted for a directory name. This is the name of your project folder. For this class, call it imd-r-training-intro. Next, you’ll select what folder to keep your project folder in. Documents/R is a good place to store all of your R projects but it’s up to you. When you are done, click on Create Project.
If you successfully started a project named imd-r-training-intro, you should see it listed at the very top right of your screen. As you start new projects, you’ll want to check that you’re working in the right one before you start coding. Take a look at the Console tab again. Notice that your current working directory is now your project folder. When you look in the Files tab of the bottom right pane, you’ll see that it also defaults to the project folder.
Let’s start coding
First we need to create a new R script file. Make sure you are working in the imd-r-training-intro project that you just created. Click on the New File icon 
in the top left corner. In the dropdown, select R Script. The source pane will appear with an untitled empty script. Go ahead and save it by clicking the Save icon (and make sure the Source on Save checkbox is deselected). Call it my_first_script.R.
The basics
We’ll start with something simple. Basic math in R is pretty straightforward and the syntax is similar to simply using a graphing calculator. Try it out! You can use the examples below or come up with your own. Even if you’re using the examples, try to actually type the code instead of copy-pasting - you’ll learn to code faster that way.
To run a single line of code in your script, place your cursor anywhere in that line and press CTRL+ENTER (or click the Run button in the top right of the script pane). To run multiple lines of code, highlight the lines you want to run and hit CTRL+ENTER or click Run.
Text following a hashtag/pound sign (#) will be ignored - use this to leave comments about what your code is doing. Commenting your code is one of the best habits you can form, and you should start now! Comments are a gift to your future self and anyone else who tries to use your code.
# By using this hashtag/pound sign, you are telling R to ignore the text afterwards. This is useful for leaving annotation of notes or instructions for yourself, or someone else using your code
# try this line to generate some basic text and become familiar with where results will appear:
print("Hello, lets do some basic math. Results of operations will appear here")## [1] "Hello, lets do some basic math. Results of operations will appear here"# one plus one
1+1## [1] 2# two times three, divided by four
(2*3)/4## [1] 1.5# basic mathematical and trigonometric functions are fairly similar to what they would be in excel
# calculate the square root of 9
sqrt(9)## [1] 3# calculate the cube root of 8 (remember that x^(1/n) gives you the nth root of x)
8^(1/3)## [1] 2# get the cosine of 180 degrees - note that trig functions in R expect angles in radians
# also note that pi is a built-in constant in R
cos(pi)## [1] -1# calculate 5 to the tenth power
5^10## [1] 9765625Notice that when you run a line of code, the code and the result appear in the console. You can also type code directly into the console, but it won’t be saved anywhere. As you get more comfortable with R, it can be helpful to use the console as a “scratchpad” for experimenting and troubleshooting. For now, it’s best to err on the side of saving your code as a script so that you don’t accidentally lose useful work.
Variables and functions
Variables
Occasionally, it’s enough to just run a line of code and display the result in the console. But typically the code we write is more complex than adding one plus one, and when we execute a line of code, we want to be able to store the result and use it later in the script. This is where variables come in. Variables allow you to assign a value (whether that’s a number, a data table, a chunk of text, or any other type of data that R can handle) to a short, human-readable name. Anywhere you put a variable in your code, R will replace it with its value when your code runs.
R uses the <- symbol for variable assignment. If you’ve used other programming languages (or remember high school algebra), you may be tempted to use = instead. It will work, but there are subtle differences between <- and =, so you should get in the habit of using <-.
# the value of 12.098 is assigned to variable 'a'
a <- 12.098
# and the value 65.3475 is assigned to variable 'b'
b <- 65.3475
# we can now perform whatever mathematical operations we want using these two variables without having to repeatedly type out the actual numbers:
a*b## [1] 790.5741(a^b)/((b+a))## [1] 7.305156e+68sqrt((a^7)/(b*2))## [1] 538.7261In the code above, we assign the variables a and b once. We can then reuse them as often as we want. This is helpful because we save ourselves some typing, reduce the chances of making a typo somewhere, and if we need to change the value of a or b, we only have to do it in one place.
Also notice that when you assign variables, you can see them listed in your Environment tab (top right pane). Remember, everything you see in the environment is just in R’s temporary memory and won’t be saved when you close out of RStudio.
All of the examples you’ve seen so far are fairly contrived for the sake of simplicity. Let’s take a look at some code that everyone here will make use of at some point: reading data from a CSV.
Functions
It’s hard to get very far in R without making use of functions. Think of a function as a black box that takes some kind of input (the argument(s)) from the user and outputs a result (the return value). Can you identify the functions in the above examples?
Using functions - reading data from CSV
The read.csv function is used to bring in a dataset from a CSV file, and it stores the data in one of the fundamental data structures in R: a data frame. read.csv() takes the file path to the CSV as input, and it outputs a data frame containing the data from the CSV.
We’ll get very familiar with data frames in this class, but for the moment just know that it’s a table of data with rows and columns. Data frames are typically organized with rows being records or observations (e.g. sampling locations, individual critters, etc.), and columns being variables that characterize those observations (e.g., species, size, date collected, x/Y coordinates). Once you have read the data in, you can take a quick look at its structure by typing the name of the variable it’s stored in.
# read in the data from tree_growth.csv and assign it as a dataframe to the variable "tree_data"
tree_data <- read.csv("data/tree_growth.csv")
# display the 'tree_data' data frame we just created
tree_data## Species DBH_in Height_ft Age_yr
## 1 arboreas madeupis 4.0 6 3
## 2 arboreas madeupis 6.0 10 5
## 3 arboreas madeupis 8.0 5 6
## 4 arboreas madeupis 11.0 7 6
## 5 arboreas madeupis 12.0 16 8
## 6 arboreas madeupis 13.0 19 8
## 7 arboreas madeupis 12.0 24 11
## 8 arboreas madeupis 20.0 22 11
## 9 arboreas madeupis 20.0 32 12
## 10 arboreas madeupis 11.0 31 13
## 11 arboreas madeupis 27.0 35 15
## 12 arboreas madeupis 25.0 33 15
## 13 arboreas madeupis 30.0 40 17
## 14 arboreas madeupis 35.0 67 18
## 15 arboreas madeupis 29.0 60 21
## 16 arboreas madeupis 37.0 65 30
## 17 arboreas madeupis 31.0 78 32
## 18 arboreas madeupis 38.0 76 36
## 19 arboreas madeupis 45.0 91 37
## 20 arboreas madeupis 55.0 83 38
## 21 arboreas madeupis 60.0 88 40
## 22 arboreas madeupis 75.0 86 43
## 23 arboreas madeupis 76.0 89 51
## 24 arboreas madeupis 80.0 75 55
## 25 arboreas madeupis 99.0 100 61
## 26 planteus falsinius 9.0 6 2
## 27 planteus falsinius 9.3 9 4
## 28 planteus falsinius 10.1 8 5
## 29 planteus falsinius 10.3 14 4
## 30 planteus falsinius 10.7 12 6
## 31 planteus falsinius 10.8 10 8
## 32 planteus falsinius 11.0 10 7
## 33 planteus falsinius 11.3 11 7
## 34 planteus falsinius 11.4 13 9
## 35 planteus falsinius 11.6 13 12
## 36 planteus falsinius 11.6 15 15
## 37 planteus falsinius 11.9 23 17
## 38 planteus falsinius 12.3 16 18
## 39 planteus falsinius 12.4 21 15
## 40 planteus falsinius 12.6 23 17
## 41 planteus falsinius 12.6 18 21
## 42 planteus falsinius 12.8 25 20
## 43 planteus falsinius 15.0 35 30Examining datasets
Compared to scrolling through a dataset in Excel, examining a dataset in R can feel unintuitive at first. Stick with it, though. Once you get used to exploring data in R, you’ll find that R will help you learn more about your dataset in a more efficient manner.
Let’s see what we can learn about the dataset as a whole. The most useful functions for this are:
head()Show only the first few lines of the dataframe in the console. The function is useful for taking a quick glance to ensure that column names were properly assigned, and to take a quick look a column names without printing the whole dataset.summary()Show a basic statistical summary of the dataset.str()Show information about the structure of the dataset, including number of rows and columns, column data types, and the first few values in each column.View()In a new tab, display a spreadsheet-like view of the full dataset with options to sort and filter. Note that you can’t change the data in this view - it is read-only.
head(tree_data) # Show the first few lines of the dataframe. Defaults to showing the first 6 rows## Species DBH_in Height_ft Age_yr
## 1 arboreas madeupis 4 6 3
## 2 arboreas madeupis 6 10 5
## 3 arboreas madeupis 8 5 6
## 4 arboreas madeupis 11 7 6
## 5 arboreas madeupis 12 16 8
## 6 arboreas madeupis 13 19 8head(tree_data, n = 10) # Show the first 10 rows## Species DBH_in Height_ft Age_yr
## 1 arboreas madeupis 4 6 3
## 2 arboreas madeupis 6 10 5
## 3 arboreas madeupis 8 5 6
## 4 arboreas madeupis 11 7 6
## 5 arboreas madeupis 12 16 8
## 6 arboreas madeupis 13 19 8
## 7 arboreas madeupis 12 24 11
## 8 arboreas madeupis 20 22 11
## 9 arboreas madeupis 20 32 12
## 10 arboreas madeupis 11 31 13summary(tree_data) # View a basic statistical summary## Species DBH_in Height_ft Age_yr
## Length:43 Min. : 4.00 Min. : 5.00 Min. : 2.00
## Class :character 1st Qu.:11.00 1st Qu.: 12.50 1st Qu.: 7.50
## Mode :character Median :12.60 Median : 23.00 Median :15.00
## Mean :24.78 Mean : 35.35 Mean :18.81
## 3rd Qu.:30.50 3rd Qu.: 62.50 3rd Qu.:25.50
## Max. :99.00 Max. :100.00 Max. :61.00str(tree_data) # Get info about the structure of the data## 'data.frame': 43 obs. of 4 variables:
## $ Species : chr "arboreas madeupis" "arboreas madeupis" "arboreas madeupis" "arboreas madeupis" ...
## $ DBH_in : num 4 6 8 11 12 13 12 20 20 11 ...
## $ Height_ft: int 6 10 5 7 16 19 24 22 32 31 ...
## $ Age_yr : int 3 5 6 6 8 8 11 11 12 13 ...View(tree_data) # Open a new tab with a spreadsheet view of the dataData structures
We’re going to take a little detour into data structures at this point. It’ll all tie back in to our tree data.
The data frame we just examined is a type of data structure. A data structure is what it sounds like: it’s a structure that holds data in an organized way. There are multiple data structures in R, including vectors, lists, arrays, matrices, data frames, and tibbles (more on this unfortunately-named data structure later). Today we’ll focus on vectors and data frames.
Vectors
Vectors are the simplest data structure in R. You can think of vectors as one column of data in an Excel spreadsheet, and the elements are each row in the column. Here are some examples of vectors:
digits <- 0:9 # Use x:y to create a sequence of integers starting at x and ending at y
digits## [1] 0 1 2 3 4 5 6 7 8 9is_odd <- rep(c(FALSE, TRUE), 5) # Use rep(x, n) to create a vector by repeating x n times
is_odd## [1] FALSE TRUE FALSE TRUE FALSE TRUE FALSE TRUE FALSE TRUEshoe_sizes <- c(7, 7.5, 8, 8.5, 9, 9.5, 10, 10.5, 11, 11.5)
shoe_sizes## [1] 7.0 7.5 8.0 8.5 9.0 9.5 10.0 10.5 11.0 11.5favorite_birds <- c("greater roadrunner", "Costa's hummingbird", "kakapo")
favorite_birds## [1] "greater roadrunner" "Costa's hummingbird" "kakapo"Note the use of c(). The c() function stands for combine, and it combines elements into a single vector. The c() function is a fairly universal way to combine multiple elements in R, and you’re going to see it over and over.
Let’s play around with vectors a little more. We can use is.vector() to test whether something is a vector. We can get the length of a vector with length(). Note that single values in R are just vectors of length one.
is.vector(digits) # Should be TRUE## [1] TRUEis.vector(favorite_birds) # Should also be TRUE## [1] TRUElength(digits) # Hopefully this is 10## [1] 10length(shoe_sizes)## [1] 10# Even single values in R are stored as vectors
length_one_chr <- "length one vector"
length_one_int <- 4
is.vector(length_one_chr)## [1] TRUEis.vector(length_one_int)## [1] TRUElength(length_one_chr)## [1] 1length(length_one_int)## [1] 1In the examples above, each vector contains a different type of data. digits contains integers, is_odd contains logical (true/false) values, favorite_birds contains text, and shoe_sizes contains decimal numbers. That’s because a given vector can only contain a single type of data. In R, there are four data types that we will typically encounter:
- character Regular text, denoted with double or single quotation marks (e.g.
"hello","3","R is my favorite programming language") - numeric Decimal numbers (e.g.
23,3.1415) - integer Integers. If you want to explicitly denote a number as an integer in R, append
Lto it or useas.integer()(e.g.5L,as.integer(30)). - logical True or false values (
TRUE,FALSE). Note thatTRUEandFALSEmust be all-uppercase.
There are two more data types, complex and raw, but you are unlikely to encounter these so we won’t cover them here.
You can use the class() function to get the data type of a vector:
class(favorite_birds)## [1] "character"class(shoe_sizes)## [1] "numeric"class(digits)## [1] "integer"class(is_odd)## [1] "logical"If you need to access a single element of a vector, you can use the syntax my_vector[x] where x is the element’s index (the number corresponding to its position in the vector). You can also use a vector of indices to extract multiple elements from the vector. Note that in R, indexing starts at 1 (i.e. my_vector[1] is the first element of my_vector). If you’ve coded in other languages, you may be used to indexing starting at 0.
second_favorite_bird <- favorite_birds[2]
second_favorite_bird## [1] "Costa's hummingbird"top_two_birds <- favorite_birds[c(1,2)]
top_two_birds## [1] "greater roadrunner" "Costa's hummingbird"Logical vectors can also be used to subset a vector. The logical vector must be the length of the vector you are subsetting.
odd_digits <- digits[is_odd]
odd_digits## [1] 1 3 5 7 9Data frames
Let’s revisit our trees data frame. We’ve explored the data frame as a whole, but it’s often useful to look at one column at a time. To do this, we’ll use the $ syntax:
tree_data$Species## [1] "arboreas madeupis" "arboreas madeupis" "arboreas madeupis"
## [4] "arboreas madeupis" "arboreas madeupis" "arboreas madeupis"
## [7] "arboreas madeupis" "arboreas madeupis" "arboreas madeupis"
## [10] "arboreas madeupis" "arboreas madeupis" "arboreas madeupis"
## [13] "arboreas madeupis" "arboreas madeupis" "arboreas madeupis"
## [16] "arboreas madeupis" "arboreas madeupis" "arboreas madeupis"
## [19] "arboreas madeupis" "arboreas madeupis" "arboreas madeupis"
## [22] "arboreas madeupis" "arboreas madeupis" "arboreas madeupis"
## [25] "arboreas madeupis" "planteus falsinius" "planteus falsinius"
## [28] "planteus falsinius" "planteus falsinius" "planteus falsinius"
## [31] "planteus falsinius" "planteus falsinius" "planteus falsinius"
## [34] "planteus falsinius" "planteus falsinius" "planteus falsinius"
## [37] "planteus falsinius" "planteus falsinius" "planteus falsinius"
## [40] "planteus falsinius" "planteus falsinius" "planteus falsinius"
## [43] "planteus falsinius"tree_data$Age_yr## [1] 3 5 6 6 8 8 11 11 12 13 15 15 17 18 21 30 32 36 37 38 40 43 51 55 61
## [26] 2 4 5 4 6 8 7 7 9 12 15 17 18 15 17 21 20 30You can also use square brackets [] to access data frame columns. You can specify columns by name or by index (integer indicating position of column). It’s almost always best to refer to columns by name when possible - it makes your code easier to read, and it prevents your code from breaking if columns get reordered.
The square bracket syntax allows you to select multiple columns at a time and to select a subset of rows.
tree_data[, "Species"] # Same as tree_data$Species## [1] "arboreas madeupis" "arboreas madeupis" "arboreas madeupis"
## [4] "arboreas madeupis" "arboreas madeupis" "arboreas madeupis"
## [7] "arboreas madeupis" "arboreas madeupis" "arboreas madeupis"
## [10] "arboreas madeupis" "arboreas madeupis" "arboreas madeupis"
## [13] "arboreas madeupis" "arboreas madeupis" "arboreas madeupis"
## [16] "arboreas madeupis" "arboreas madeupis" "arboreas madeupis"
## [19] "arboreas madeupis" "arboreas madeupis" "arboreas madeupis"
## [22] "arboreas madeupis" "arboreas madeupis" "arboreas madeupis"
## [25] "arboreas madeupis" "planteus falsinius" "planteus falsinius"
## [28] "planteus falsinius" "planteus falsinius" "planteus falsinius"
## [31] "planteus falsinius" "planteus falsinius" "planteus falsinius"
## [34] "planteus falsinius" "planteus falsinius" "planteus falsinius"
## [37] "planteus falsinius" "planteus falsinius" "planteus falsinius"
## [40] "planteus falsinius" "planteus falsinius" "planteus falsinius"
## [43] "planteus falsinius"tree_data[, c("Species", "Age_yr")] # Data frame with only Species and Age_yr columns## Species Age_yr
## 1 arboreas madeupis 3
## 2 arboreas madeupis 5
## 3 arboreas madeupis 6
## 4 arboreas madeupis 6
## 5 arboreas madeupis 8
## 6 arboreas madeupis 8
## 7 arboreas madeupis 11
## 8 arboreas madeupis 11
## 9 arboreas madeupis 12
## 10 arboreas madeupis 13
## 11 arboreas madeupis 15
## 12 arboreas madeupis 15
## 13 arboreas madeupis 17
## 14 arboreas madeupis 18
## 15 arboreas madeupis 21
## 16 arboreas madeupis 30
## 17 arboreas madeupis 32
## 18 arboreas madeupis 36
## 19 arboreas madeupis 37
## 20 arboreas madeupis 38
## 21 arboreas madeupis 40
## 22 arboreas madeupis 43
## 23 arboreas madeupis 51
## 24 arboreas madeupis 55
## 25 arboreas madeupis 61
## 26 planteus falsinius 2
## 27 planteus falsinius 4
## 28 planteus falsinius 5
## 29 planteus falsinius 4
## 30 planteus falsinius 6
## 31 planteus falsinius 8
## 32 planteus falsinius 7
## 33 planteus falsinius 7
## 34 planteus falsinius 9
## 35 planteus falsinius 12
## 36 planteus falsinius 15
## 37 planteus falsinius 17
## 38 planteus falsinius 18
## 39 planteus falsinius 15
## 40 planteus falsinius 17
## 41 planteus falsinius 21
## 42 planteus falsinius 20
## 43 planteus falsinius 30tree_data[1:5, "Species"] # First five elements of Species column## [1] "arboreas madeupis" "arboreas madeupis" "arboreas madeupis"
## [4] "arboreas madeupis" "arboreas madeupis"Now that we’ve learned about vectors, these data frame columns might look familiar. A data frame is just a collection of vectors of the same length, where each vector is a column.
is.vector(tree_data$Species)## [1] TRUEclass(tree_data$Species)## [1] "character"is.vector(tree_data$Age_yr)## [1] TRUEclass(tree_data$Age_yr)## [1] "integer"str(tree_data) # Check out the abbreviations next to the column names: chr, num, and int. These are the data types of the columns.## 'data.frame': 43 obs. of 4 variables:
## $ Species : chr "arboreas madeupis" "arboreas madeupis" "arboreas madeupis" "arboreas madeupis" ...
## $ DBH_in : num 4 6 8 11 12 13 12 20 20 11 ...
## $ Height_ft: int 6 10 5 7 16 19 24 22 32 31 ...
## $ Age_yr : int 3 5 6 6 8 8 11 11 12 13 ...Working with columns
R is great at working with vectors because they are the fundamental data structure of the language. Since data frame columns are vectors, they are typically easy to operate on.
Basic statistical summary functions aren’t too different from what you’d do in Excel.
#to calculate the mean of the 'age' column in the original dataframe:
mean(tree_data$Age_yr)## [1] 18.81395#or to calculate the mean of the DBH vector we created:
dbh <- tree_data$DBH_in
mean(dbh)## [1] 24.78372#like-wise for the median, standard deviation, and minimum:
median(tree_data$Age_yr)## [1] 15sd(tree_data$Age_yr)## [1] 15.02261min(tree_data$Age_yr)## [1] 2Another useful way to summarize a column is to get the unique values it contains. This is typically most useful with character columns (e.g. species) but works with any data type.
unique(tree_data$Species) # in this line, we are printing all of the unique species names in the dataset (in this case 2).## [1] "arboreas madeupis" "planteus falsinius"species_names <- unique(tree_data$Species) # assign unique species names to a variableSummarizing data is useful but we also want to be able to modify it. So let’s say we want to convert the ‘Height_ft’ column in tree_data from feet to inches. We’ll start by assigning a conversion factor to the variable in_per_ft.
in_per_ft <- 12 # Since there are 12 inches in a footWe can then use this conversion factor to convert the ‘Height_ft’ column to inches, and place this converted value in a new column we will create in the data frame. Notice below that if you assign a vector to a column that doesn’t exist yet, R will automatically create that column. We’ll use the head() function to verify that the Height_in column was created correctly.
# here we are specifying the new column on the left side of the '<-', and telling R what we want put into it on the right side of the '<-'.
tree_data$Height_in <- tree_data$Height_ft * in_per_ft
# we can now use the 'head function to check our work, and make sure the proper conversion was carried out:
head(tree_data)## Species DBH_in Height_ft Age_yr Height_in
## 1 arboreas madeupis 4 6 3 72
## 2 arboreas madeupis 6 10 5 120
## 3 arboreas madeupis 8 5 6 60
## 4 arboreas madeupis 11 7 6 84
## 5 arboreas madeupis 12 16 8 192
## 6 arboreas madeupis 13 19 8 228Basic data visualization
It’s often easier to get a sense of our data by visualizing it, so let’s explore the basic tools for taking a first glance at our data.
Histogram
We’ll start with a histogram. The code below generates a basic histogram plot of a specific column in the dataframe (remember from earlier that this is denoted by the $ sign) using the hist() function. We’ll start by creating a histogram for the ‘Age_yr’ column.
hist(x = tree_data$Age_yr)
Looking at the histogram, we can see that the age of most trees in the dataset fall in the 0-20 year old age class.
Scatter plot
The data frame also includes information on DBH, and height of the trees. What if we want to quickly take a look at how these two variables relate to one another? To get a basic plot, you can use the plot() function. A simple call to plot() takes two arguments: x (a vector of x coordinates) and y (a vector of y coordinates). You can learn more about plot() by typing ?plot at the console.
# this should generate an ugly but effective scatterplot of tree height vs age. It would appear that older trees are taller!
plot(x = tree_data$Age_yr, y = tree_data$Height_ft)
# let's try the same, but with DBH as the dependent variable on the y axis:
# again, we should see a plot below showing that, even in a make-believe forest, older trees of both species tend to be thicker!
plot(x = tree_data$Age_yr, y = tree_data$DBH_in)
Getting Help
Getting Help
There are a number of options to get help with R. If you’re trying to figure out how to use a function, you can type ?function_name. For example ?plot will show the R documentation for that function in the Help panel.
?plot
?dplyr::filterIf you have a question, online resources include Stackexchange and Stackoverflow are extremely helpful. Google searches are a great first step. It helps to include “in R” in every search related to R code.
Don’t hesitate to reach out to colleagues for help as well! If you are stuck on something and the answers on Google are more confusing than helpful, don’t be afraid to ask a human. Every experienced R programmer was a beginner once, so chances are they’ve encountered the same problem as you at some point. There is an R-focused data science community of practice for I&M folks, which anyone working in R (regardless of experience!) is invited and encouraged to join. The Teams join code will be shared with course participants.
Day 2: Data Wrangling
Intro and Packages
Data Wrangling Intro
Goals of this Day
- Learn to load and understand how packages work
- Understand what data wrangling is and some ways to do it.
- Learn how to clean your import data so that R understands what it is
- Get data in the format you need for analysis using the tidyverse.
- Combine data from different data.frames so that it is in the correct format
Disclaimer: While R provides lots of ways to check and clean up data, these tools are no replacement for sound data management that prevents data issues from happening in the first place.
Working with Packages: Installing a Package
What is a package? A package is a group of functions that someone has put together. They come with documentation that tells you about what each function does.
First, load package
- Go to packages tab in lower right quadrant
- Click on “install”
- Type name of package ‘janitor’
- Click “install”
- Code will pop up in your console
install.packages("janitor")Calling a Package
To use a package, call it by using the function “library”
library('janitor')This is called “loading a package” or “calling a package”.
Using Package
Look at the vignette
- Go to packages window of lower right quadrant
- Click on ‘janitor’ name in list of packages
- Vignette will show up in the ‘help’ tab
take time to click a function from the vignette list walk through how to get helpful info
Did you get “error” text?
Text pops up in the console after you run code You must decide whether it means you have to change something or not
library('janitor')##
## Attaching package: 'janitor'## The following objects are masked from 'package:stats':
##
## chisq.test, fisher.testChallenge: Install a Package
Challenge: install + call the package ‘tidyverse’
Solution
Solution: install + call the package ‘tidyverse’
install.packages("tidyverse")
library('tidyverse')Data Cleaing
Data Cleaning
Data wrangling is all the steps you have to take to prepare your data for analysis or plotting. Data munging or data cleaning are other terms used to describe this process.Material presented today can be found in Hadley Wickham’s book “R for Data Science”, which is free and online. https://r4ds.had.co.nz/wrangle-intro.html
If you’re starting from this tab of the lesson, you’ll need to load two packages before proceeding.
library('janitor')
library('tidyverse')Major Components to Data Wrangling
- Step 1: Import (covered yesterday)
- Step 2+3: Tidy & Transform (get data into a useful format)
Column Names
Sometimes you get data where the column names are messy. For an example of data that needs tidying before exploration and analyses, we’ll use fake data. Copy + paste the code chunk below into your script, and run it in the console to generate a dataframe called “messy_data”.
Here, we’ll clarify that we are using the word “data.frame” or “dataframe” to describe the data. The word “tibble” also shows up. We’ll use these terms interchangeably.
messy_data <- tibble(TrailName = c('Bayside', 'Coastal', 'Tidepools',
'Monument', 'Tidepools', NA),
Visitor_count = c(200, 300, 400, 500, 400, NA),
`empty row` = NA)This data.frame captures the hypothetical number of visitors across four trails at a park. There are a few qualities that make it difficult to interpret. Seeing this data, what would you change to make it most consistent?
messy_data## # A tibble: 6 x 3
## TrailName Visitor_count `empty row`
## <chr> <dbl> <lgl>
## 1 Bayside 200 NA
## 2 Coastal 300 NA
## 3 Tidepools 400 NA
## 4 Monument 500 NA
## 5 Tidepools 400 NA
## 6 <NA> NA NAThere are four things I suggest fixing with these data:
- The column names are inconsistently formatted
- There is an empty row
- There is an empty column
- There is a duplicated row
The Janitor package has ways to fix all of these data situations and aid in data cleaning. Let’s proceed with changing the column names.
The clean_names() function in the janitor package will consistently format column names.
clean_data <- clean_names(messy_data)
clean_data## # A tibble: 6 x 3
## trail_name visitor_count empty_row
## <chr> <dbl> <lgl>
## 1 Bayside 200 NA
## 2 Coastal 300 NA
## 3 Tidepools 400 NA
## 4 Monument 500 NA
## 5 Tidepools 400 NA
## 6 <NA> NA NAEmpty Rows and Columns
Now, we’ll deal with the empty row and empty column. To remove empty rows or columns, we can use the remove_empty() function.
The arguments for the function are as follows:
remove_empty(
# first argument is the data
dat,
# the next specifies if you want to get rid of empty rows, columns, or both
which = c("rows", "cols"),
# this last one asks whether (TRUE) or not (FALSE) you want R to tell you which rows or columns were removed
quiet = TRUE)Let’s give it a try
Removing empty rows
clean_data2 <- remove_empty(clean_data, which = c('rows'), quiet = FALSE)
clean_data2## # A tibble: 5 x 3
## trail_name visitor_count empty_row
## <chr> <dbl> <lgl>
## 1 Bayside 200 NA
## 2 Coastal 300 NA
## 3 Tidepools 400 NA
## 4 Monument 500 NA
## 5 Tidepools 400 NARemoving empty rows and columns
clean_data2 <- remove_empty(clean_data, which = c('rows', 'cols'), quiet = TRUE)
clean_data2## # A tibble: 5 x 2
## trail_name visitor_count
## <chr> <dbl>
## 1 Bayside 200
## 2 Coastal 300
## 3 Tidepools 400
## 4 Monument 500
## 5 Tidepools 400The Pipe - Sequences of Functions
What happens if you’d like to both clean_names() and remove_empty() on the same data?
Solution 1 - run as a sequence
clean_data <- clean_names(messy_data)
clean_data <- remove_empty(clean_data, which = c('rows', 'cols'), quiet = TRUE)
clean_data## # A tibble: 5 x 2
## trail_name visitor_count
## <chr> <dbl>
## 1 Bayside 200
## 2 Coastal 300
## 3 Tidepools 400
## 4 Monument 500
## 5 Tidepools 400Solution 2 - use the pipe operator
Pipe operators (%>% or |>) can be read like the word “then”. For example, clean_names() %>% remove_empty() can be read as “clean names, then remove empty”. They allow for a more clean workflow.
clean_data <- clean_names(messy_data) %>%
remove_empty(which = c('rows', 'cols'), quiet = TRUE)
clean_data## # A tibble: 5 x 2
## trail_name visitor_count
## <chr> <dbl>
## 1 Bayside 200
## 2 Coastal 300
## 3 Tidepools 400
## 4 Monument 500
## 5 Tidepools 400# note that replacing %>% with |> also worksNote that there is no data argument in the remove_empty() function. This is because, when using a pipe, the data argument is inherited from the prior function.
Removing Duplicate Entries
A final function we’ll discuss is distinct() from the dplyr package, which is park of tidyverse. It removes duplicate rows from data.
clean_data <- clean_names(messy_data) %>%
remove_empty(which = c('rows', 'cols'), quiet = TRUE) %>%
distinct()
clean_data## # A tibble: 4 x 2
## trail_name visitor_count
## <chr> <dbl>
## 1 Bayside 200
## 2 Coastal 300
## 3 Tidepools 400
## 4 Monument 500A word of caution - sometimes you may not want to remove duplicate rows, or may not want to remove empty rows, because they are important to the data. You know your data best. Use functions responsibly.
Challenge
The code chunk below creates messy data, called “messy_data2”. Use functions from the Janitor package to tidy the data, and name the result “clean_data2”. See if you can use the pipe operator to streamline the script.
messy_data2 <- tibble(`Park Code` = c(NA, NA, 'CABR', 'CHIS', 'SAMO', 'SAMO'),
visitor_count = c(NA, NA, 400, 500, 400, 400),
`empty row` = NA)Solution
clean_data2 <- clean_names(messy_data2) %>%
remove_empty(which = c('rows', 'cols'), quiet = TRUE) %>%
distinct()Dplyr functions
Data Wrangling
In this section, we’ll introduce functions in the dplyr package. This package is one of the packages that make up the tidyverse, so you won’t have to install and call another package. The dplyr package is a group of functions that help to manipulate and summarize data. https://www.rdocumentation.org/packages/dplyr/versions/0.7.8
You can think of these functions like verbs that accomplish tasks, and are strung together by the pipe operator (%>% or |>).
We’ll cover the following functions:
- select() - selects or removes columns
- rename() - renames columns
- arrange() - arranges a data.frame by column values
- filter() - filters data for observations that meet specific criteria
- mutate() - adds new column(s) to a data.frame
- pull() %>% table() - isolates and summarizes by a column
- group_by() %>% summarize() - summary calculations by groups/factor levels and returns data for each group
- group_by() %>% mutate() - summary calculations by groups/factor levels and returns data for each observation
The data we’ll use for modules 3-5 is from the Tidy Tuesday community on GitHub. Each week on Tuesday, the community posts data for exploratory data analysis and visualization. You can add data from this page to your RStudio environment by copying + pasting the “Get the data!” code chunk into your script and running it in the console.
https://github.com/rfordatascience/tidytuesday/tree/master/data/2019/2019-09-17
The data we’ll use is in csv format, so we’ll need to call the readr package, as well as the tidyverse and janitor for this lesson (if you’ve already called tidyverse and janitor, don’t run those lines).
library('readr')
library('janitor')
library('tidyverse')Next, we’ll load the data from the R for Data Science Tidy Tuesday GitHub page. There are three lines which will load three data.frames.
- park_visits - which contains data on visits to US National Parks
- state_pop - which contains information on populations of US states
- gas_price - which contains information on gas prices over space and time
For now, we’ll focus on the park_visits data.frame, and will circle back to the other two later. There are descriptions of each column for each data.frame on the Tidy Tuesday website linked above if you’d like more information.
park_visits <- read_csv("https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2019/2019-09-17/national_parks.csv")
state_pop <- read_csv("https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2019/2019-09-17/state_pop.csv")
gas_price <- read_csv("https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2019/2019-09-17/gas_price.csv")Working with dplyr
Dplyr is a package that is part of the tidyverse, which is used with data.frames. It makes tasks like filtering and summarizing data easier. Many of the functions within the package are “verb-like” so they are easy to understand and remember.
dplyr - select() and rename()
select() - keeps the columns you indicate and drops the rest. Since select() puts the columns in the order you type them, you can also use select() to reorder columns.
rename() - keeps all the columns, but renames the ones you indicate
Since we won’t be using geographical information, other than state, in our module, we’ll select some columns from the park_visits data to use today.
Note that, the name to the left of the arrow is the resulting data.frame, the information to the right that’s connected by pipes are the steps to get to the result. I normally name the result something different than the original data.frame so that the original data aren’t affected. However, since we won’t need the geographical information later on, I will name this one “park_visits”.
# let's see the column names in park_visits
colnames(park_visits)## [1] "year" "gnis_id" "geometry"
## [4] "metadata" "number_of_records" "parkname"
## [7] "region" "state" "unit_code"
## [10] "unit_name" "unit_type" "visitors"# we can also look at the first few rows by using head()
head(park_visits)## # A tibble: 6 x 12
## year gnis_id geometry metadata number_of_records parkname region state
## <chr> <dbl> <chr> <chr> <dbl> <chr> <chr> <chr>
## 1 1904 1163670 POLYGON <NA> 1 Crater Lake PW OR
## 2 1941 1531834 MULTIPOLYGON <NA> 1 Lake Roose~ PW WA
## 3 1961 2055170 MULTIPOLYGON <NA> 1 Lewis and ~ PW WA
## 4 1935 1530459 MULTIPOLYGON <NA> 1 Olympic PW WA
## 5 1982 277263 POLYGON <NA> 1 Santa Moni~ PW CA
## 6 1919 578853 MULTIPOLYGON <NA> 1 <NA> NE ME
## # ... with 4 more variables: unit_code <chr>, unit_name <chr>, unit_type <chr>,
## # visitors <dbl># you can either select for the columns you want
park_visits1 <- park_visits %>% select(year, parkname, unit_code, state, visitors)
# or against the columns you don't want, by using the minus sign before the column name
park_visits2 <- park_visits %>% select(-gnis_id, -geometry, -metadata, -number_of_records, -region, -unit_name, -unit_type)
# you can check if you have the same result by using colnames()
colnames(park_visits1)## [1] "year" "parkname" "unit_code" "state" "visitors"colnames(park_visits2)## [1] "year" "parkname" "state" "unit_code" "visitors"# to clear up our environment, we'll get rid of park_visits1 and 2 using remove()
remove(park_visits1, park_visits2)
# since we won't need all of the columns in park_visits, we'll change that file (note that park_visits is on both the left and right sides of the pipe)
park_visits <- park_visits %>% select(year, parkname, unit_code, state, visitors)Taking a look at the park_visits data, there are some column name inconsistencies. For example, some are in snake case (unit_code), but some are not (parkname). To rename a column, use rename(). The structure of this argument is as follows:
rename(dataset, new_name = old_name) OR dataset %>% rename(new_name = old_name)
# rename parkname column to make it snake case
park_visits <- park_visits %>% rename(park_name = parkname)
# we can check if we've done it correctly by looking at the column names
colnames(park_visits)## [1] "year" "park_name" "unit_code" "state" "visitors"dplyr- filter()
filter() - filters out rows from a data.frame based on values in columns can combine criteria using & (and), | (or), ! (not)
If you’d like to get visitation data for just parks in the state of Washington (WA), you can use the filter function. The structure is as follows:
wa_park_visits <- filter(
# name of data.frame
park_visits,
# condition
state == 'WA'
)
# head() shows the first few rows of data
head(wa_park_visits)
# you can use the unique() to check you've done the filtering correctly
unique(wa_park_visits$state)Now that you’ve got parks in the state of Washington, perhaps you’d like to know which parks in WA had over 5 million total visitors. This will require more intensive filtering. We’ll have to filter for parks where state == ‘WA’ AND year == ‘Total’ AND visitation >= 5,000,000.
wa_park_visits <- park_visits %>%
# filter for parks in WA
filter(state == 'WA' &
# filter for total visitation
year == 'Total' &
# filter for over 5 million visitors
visitors > 5000000)
# head() shows the first few rows of data
head(wa_park_visits)You also might be wondering why we use “==” instead of “=”. In R, the single equal sign “=” is used to assign something. For example, if you want to make some variable equal to a number, you could tell R “x = 5”, and it would remember that x is equal to 5. To evaluate something as TRUE/FALSE, you use “==”. If you type “5 == 3” into the console, it will return the answer “FALSE”. The filtering argument uses this TRUE/FALSE evaluation - filtering for rows where the stated condition - such as year == ’Total” - is TRUE.
dplyr- challenge 1
Time for the first challenge. Using the park_visits data, create a new data.frame containing total visitation for a state of your choice. Note that some “states” aren’t states. For example, VI refers to the Virgin Islands. How can you check your work?
dplyr- solution 1
pr_park_visits <- park_visits %>%
# filter for parks in Puerto Rico
filter(state == 'PR' &
# filter for total visitation
year == 'Total')
# head() shows the first few rows of data
head(pr_park_visits)
# San Juan National Historic Site (Puerto Rico) had > 55 million visitors!
# remove result from environment
remove(pr_park_visits)dplyr - arrange()
arrange() reorders rows. You can think of this like “sort” in Excel. While this is generally not needed in R, it can be helpful for exploring your data, and when working with time series data and lagged values.
As an example, we’ll re-visit total park visitation in WA, sorting by the number of visitors.
wa_park_visits <- park_visits %>%
# filter for parks in WA
filter(state == 'WA' &
# filter for total visitation
year == 'Total' &
# filter for over 5 million visitors
visitors > 5000000) %>%
# arrange by visitation - default is ascending order
arrange(visitors)
# look at result (df is short so head() not needed)
wa_park_visits## # A tibble: 9 x 5
## year park_name unit_code state visitors
## <chr> <chr> <chr> <chr> <dbl>
## 1 Total Lewis and Clark LEWI WA 9515314
## 2 Total <NA> NEPE WA 9566104
## 3 Total Ross Lake ROLA WA 11335924
## 4 Total North Cascades NOCA WA 14809679
## 5 Total Fort Vancouver FOVA WA 20238294
## 6 Total San Juan Island SAJH WA 23066952
## 7 Total Lake Roosevelt LARO WA 62247752
## 8 Total Mount Rainier MORA WA 94488801
## 9 Total Olympic OLYM WA 160737962# to arrange by descending order, we can use desc()
wa_park_visits2 <- wa_park_visits %>%
arrange(desc(visitors))
# look at result
wa_park_visits2## # A tibble: 9 x 5
## year park_name unit_code state visitors
## <chr> <chr> <chr> <chr> <dbl>
## 1 Total Olympic OLYM WA 160737962
## 2 Total Mount Rainier MORA WA 94488801
## 3 Total Lake Roosevelt LARO WA 62247752
## 4 Total San Juan Island SAJH WA 23066952
## 5 Total Fort Vancouver FOVA WA 20238294
## 6 Total North Cascades NOCA WA 14809679
## 7 Total Ross Lake ROLA WA 11335924
## 8 Total <NA> NEPE WA 9566104
## 9 Total Lewis and Clark LEWI WA 9515314# remove one result
remove(wa_park_visits2)dplyr - pull()
pull() takes data from a column and returns in in vector form
pull() %>% table() is a handy way to see a summary of categorical data (e.g. how may observations for each state there are in your data)
# two ways to get the number of observations per state
# option 1 - use pull() to isolate a column and table() to get # of observations
state_n <- park_visits %>%
# pull the state column only
pull(state) %>%
# get # of observations by state in table output
table()
# option 2 - use group_by() to isolate column(s) and tally() to get # of obs
state_n2 <- park_visits %>%
# group by state
group_by(state) %>%
# get tally of observations by state in output
tally()
# remove created objects
remove(state_n, state_n2)dplyr - mutate()
mutate() makes new columns by doing calculations on old columns
As an example, let’s calculate the total visitation for parks with over 5 million total visitors in WA in millions, by dividing the number of total visitors by 1,000,000.
wa_park_visits_millions <- wa_park_visits %>%
mutate(visitors_mil = visitors/1000000)
wa_park_visits_millions## # A tibble: 9 x 6
## year park_name unit_code state visitors visitors_mil
## <chr> <chr> <chr> <chr> <dbl> <dbl>
## 1 Total Lewis and Clark LEWI WA 9515314 9.52
## 2 Total <NA> NEPE WA 9566104 9.57
## 3 Total Ross Lake ROLA WA 11335924 11.3
## 4 Total North Cascades NOCA WA 14809679 14.8
## 5 Total Fort Vancouver FOVA WA 20238294 20.2
## 6 Total San Juan Island SAJH WA 23066952 23.1
## 7 Total Lake Roosevelt LARO WA 62247752 62.2
## 8 Total Mount Rainier MORA WA 94488801 94.5
## 9 Total Olympic OLYM WA 160737962 161.dplyr - group_by() and mutate()
group_by() indicates that data is in groups, so you can do calculations separately for each group
group_by() %>% mutate() will add a column and return the same number of rows as the original data
This is helpful when you want to get group summary information without sacrificing the original values. For example, if we want to know what percent of the total visitation of a park occurs in each year, and add this to the original data.
park_visits2 <- park_visits %>%
# remove the rows with Totals per park
filter(year != "Total") %>%
# do calculations by park
group_by(park_name) %>%
# add visit_percent column
# visits fore each year divided by the sum of total visits for each park
# the na.rm = TRUE means that NA values are not included in calculations
# round(2) is rounds result to 2 decimal places for easier reading
mutate(visit_percent = (100*visitors/sum(visitors, na.rm = TRUE)) %>%
round(2))
# take a look at the result
head(park_visits2)## # A tibble: 6 x 6
## # Groups: park_name [6]
## year park_name unit_code state visitors visit_percent
## <chr> <chr> <chr> <chr> <dbl> <dbl>
## 1 1904 Crater Lake CRLA OR 1500 0
## 2 1941 Lake Roosevelt LARO WA 0 0
## 3 1961 Lewis and Clark LEWI WA 69000 0.73
## 4 1935 Olympic OLYM WA 2200 0
## 5 1982 Santa Monica Mountains SAMO CA 468144 2.66
## 6 1919 <NA> ACAD ME 64000 0# remove from environment
remove(park_visits2)dplyr - group_by() and summarize()
group_by() and summarize() (or summarise()) also work with grouped data. Here you calculate a summary for each group. Output data frame will have the columns that define the group and the summary data, will have one row for each group.
I find this combination most helpful for calculating means and standard deviations for data. You can also find minimum and maximum values. We’ll filter for total visitation numbers for each park, then calculate a visitation summary for each state.
state_summary <- park_visits %>%
# filter for total visitation numbers (to avoid double counts)
filter(year == 'Total') %>%
# do calculations by state
group_by(state) %>%
# calculate summaries
summarize(
# mean number of total visitors across parks in each state
mean_visit = mean(visitors, na.rm = TRUE),
# sd of total visitors across parks in each state
sd_visit = sd(visitors, na.rm = TRUE),
# min total visitors across parks in each state
min_visit = min(visitors, na.rm = TRUE),
# max total visitors across parks in each state
max_visit = max(visitors, na.rm = TRUE),
# get number of park units w data in each state
n_parks = n()
)
# take a look at the result
head(state_summary)## # A tibble: 6 x 6
## state mean_visit sd_visit min_visit max_visit n_parks
## <chr> <dbl> <dbl> <dbl> <dbl> <int>
## 1 AK 4661539 6737110. 10490 19476522 22
## 2 AL 2998848. 2162876. 365883 5487784 5
## 3 AR 20647567. 34237894. 60563 92028933 7
## 4 AS 113694 NA 113694 113694 1
## 5 AZ 23655713. 46310344. 348811 205486894 19
## 6 CA 61322519. 120602434. 6349 611031225 26# if you want to continue analyzing this summarized data, you must first ungroup()
# for example:
state_summary <- state_summary %>%
ungroup() %>%
mutate(dif = max_visit - min_visit)
# remove from environment
remove(state_summary)dplyr - challenge 2
For this second challenge, you’ll start with the park_visits data. Calculate the average visitation divided by the number of parks for each state. The final data.frame should be arranged in descending order by a numeric column of your choice.
dplyr - solution 2
challenge2 <- park_visits %>%
# filter for total visitation numbers (to avoid double counts)
filter(year == 'Total') %>%
# do calculations by state
group_by(state) %>%
# calculate summaries
summarize(
# mean number of total visitors across parks in each state
mean_visit = mean(visitors, na.rm = TRUE),
# get number of park units w data in each state
n_parks = n()
) %>%
# ungroup to do another calculation
ungroup() %>%
# calculate visit/n
mutate(visit_per_park = mean_visit/n_parks) %>%
# select relevant columns
select(state, visit_per_park) %>%
# arrange in descending order
arrange(desc(visit_per_park))
# take a look at the result
head(challenge2)## # A tibble: 6 x 2
## state visit_per_park
## <chr> <dbl>
## 1 NV 103867788.
## 2 PR 55928192
## 3 ME 40890226.
## 4 NC 19591218.
## 5 MS 19160525
## 6 OK 18206483.# remove from environment
remove(challenge2)Pivoting Tables
Pivoting Tables
This section will go over the use of functions to change data from wide to long format and vice-versa.
We will:
- Explain what these formats are and why anyone cares
pivot_longer()pivot_Wider()
If you’re starting with this section, run the following code chunk to get caught up.
# call packages
library("readr")
library("janitor")
library("tidyverse")
# get data
park_visits <- read_csv("https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2019/2019-09-17/national_parks.csv")
state_pop <- read_csv("https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2019/2019-09-17/state_pop.csv")
gas_price <- read_csv("https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2019/2019-09-17/gas_price.csv")
# cut some columns from park_visits data
park_visits <- park_visits %>% select(year, parkname, unit_code, state, visitors)Why do we care?
Depending on your goal, you may want to have your data formatted in multiple ways. For example, it is much easier for humans to read wide format data, while most R functions rely on long format data. Here is a quick demonstration using park visits data:
# we'll first take a subset of the park visits data
medn_visits <- park_visits %>%
# get parks in MEDN and no totals, years 2010-2015
filter(unit_code %in% c("CABR", "CHIS", "SAMO") & year != "Total") %>%
# make year a number, since no more text
mutate(year = as.numeric(year)) %>%
# get years 2010:2015
filter(year >= 2010 & year <= 2015) %>%
# arrange in ascending order
arrange(year) %>%
# select relevant columns
select(unit_code, year, visitors)
# if we take a look at the data, it is in long format, each observation is a row
# this isn't the most human-friendly way to look at data, but it is really helpful if we want to use dplyr group_by() functions like mutate() and summarize()
medn_visits## # A tibble: 18 x 3
## unit_code year visitors
## <chr> <dbl> <dbl>
## 1 SAMO 2010 568371
## 2 CABR 2010 763140
## 3 CHIS 2010 277515
## 4 SAMO 2011 609636
## 5 CABR 2011 813351
## 6 CHIS 2011 242756
## 7 SAMO 2012 649471
## 8 CABR 2012 877951
## 9 CHIS 2012 249594
## 10 SAMO 2013 633054
## 11 CABR 2013 856128
## 12 CHIS 2013 212029
## 13 SAMO 2014 694714
## 14 CABR 2014 893434
## 15 CHIS 2014 342161
## 16 SAMO 2015 797126
## 17 CABR 2015 981825
## 18 CHIS 2015 324816# conversely, we can pivot longer to make the data wide format
# it is more human-friendly but you can't use dplyr functions on it easily
medn_visits_wide <- pivot_wider(medn_visits, names_from = year, values_from = visitors)
medn_visits_wide## # A tibble: 3 x 7
## unit_code `2010` `2011` `2012` `2013` `2014` `2015`
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 SAMO 568371 609636 649471 633054 694714 797126
## 2 CABR 763140 813351 877951 856128 893434 981825
## 3 CHIS 277515 242756 249594 212029 342161 324816There are two functions we’ll discuss today:
pivot_wider(): converts from long –> widepivot_longer(): converts from wide –> long
pivot_wider()
Full disclosure, I almost always have to google the arguments that go into the pivot_ functions. Here’s the website for pivot_wider() https://tidyr.tidyverse.org/reference/pivot_wider.html
Here’s a short example - pivoting the MEDN visits from long to wide, with more detail than the prior code chunk.
# check out the original data (medn_visits) in long format
medn_visits## # A tibble: 18 x 3
## unit_code year visitors
## <chr> <dbl> <dbl>
## 1 SAMO 2010 568371
## 2 CABR 2010 763140
## 3 CHIS 2010 277515
## 4 SAMO 2011 609636
## 5 CABR 2011 813351
## 6 CHIS 2011 242756
## 7 SAMO 2012 649471
## 8 CABR 2012 877951
## 9 CHIS 2012 249594
## 10 SAMO 2013 633054
## 11 CABR 2013 856128
## 12 CHIS 2013 212029
## 13 SAMO 2014 694714
## 14 CABR 2014 893434
## 15 CHIS 2014 342161
## 16 SAMO 2015 797126
## 17 CABR 2015 981825
## 18 CHIS 2015 324816# if we want to make each row a park and each column a year, we can pivot wider
medn_visits_wide <- medn_visits %>%
pivot_wider(
# names from is the thing you want to be the columns
names_from = year,
# values from is what you want to fill the columns with
values_from = visitors
)
# check out the result
medn_visits_wide## # A tibble: 3 x 7
## unit_code `2010` `2011` `2012` `2013` `2014` `2015`
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 SAMO 568371 609636 649471 633054 694714 797126
## 2 CABR 763140 813351 877951 856128 893434 981825
## 3 CHIS 277515 242756 249594 212029 342161 324816# we can also make park units the columns and fill it with visitor data for each year
medn_visits_wide <- medn_visits %>%
pivot_wider(names_from = unit_code, values_from = visitors)
medn_visits_wide## # A tibble: 6 x 4
## year SAMO CABR CHIS
## <dbl> <dbl> <dbl> <dbl>
## 1 2010 568371 763140 277515
## 2 2011 609636 813351 242756
## 3 2012 649471 877951 249594
## 4 2013 633054 856128 212029
## 5 2014 694714 893434 342161
## 6 2015 797126 981825 324816pivot_longer()
Now that we have our data in wide format, it is time to revert back to pivoting longer. Here’s the reference material for pivot_longer() https://tidyr.tidyverse.org/reference/pivot_longer.html.
# check out the original data (medn_visits_wide) in wide format
medn_visits_wide## # A tibble: 6 x 4
## year SAMO CABR CHIS
## <dbl> <dbl> <dbl> <dbl>
## 1 2010 568371 763140 277515
## 2 2011 609636 813351 242756
## 3 2012 649471 877951 249594
## 4 2013 633054 856128 212029
## 5 2014 694714 893434 342161
## 6 2015 797126 981825 324816# if we want to make each row an observation, with an associated park and year, we'll pivot_longer
medn_visits_long <- medn_visits_wide %>%
pivot_longer(
# first argument is the columns you'd like to transform
SAMO:CHIS,
# next is what you'd like to name the former name cols
names_to = "unit_code",
# last is what you'd like to name the former values cols
values_to = "visitors"
)
# check out the result
medn_visits_long## # A tibble: 18 x 3
## year unit_code visitors
## <dbl> <chr> <dbl>
## 1 2010 SAMO 568371
## 2 2010 CABR 763140
## 3 2010 CHIS 277515
## 4 2011 SAMO 609636
## 5 2011 CABR 813351
## 6 2011 CHIS 242756
## 7 2012 SAMO 649471
## 8 2012 CABR 877951
## 9 2012 CHIS 249594
## 10 2013 SAMO 633054
## 11 2013 CABR 856128
## 12 2013 CHIS 212029
## 13 2014 SAMO 694714
## 14 2014 CABR 893434
## 15 2014 CHIS 342161
## 16 2015 SAMO 797126
## 17 2015 CABR 981825
## 18 2015 CHIS 324816# this should look familiar - like the original medn_visits datapivot challenge
This is a more involved challenge that will draw on your knowledge of both dplyr (from Module 3) and pivoting.
Using the state_pop data, accomplish the following tasks:
- Get data for 3 states of your choice for the years 2010-2015 (note - there is no ‘total’ values in the year column)
- Create a
data.framein long format - Pivot wider to make a human-friendly table with states as rows and years as columns
- Pivot longer, returning to the original layout
Hint: if your column names are numbers, like years, you must refer to them using back-ticks. A column called 2020 would be `2020`
pivot solution
# get 3 states of choice (FL, CA, AK) for years 2010-2015
solution_original <- state_pop %>%
filter(state %in% c("FL", "AK", "CA") &
year >= 2010 & year <= 2015)
# pivot wider
solution_wide <- solution_original %>%
pivot_wider(
names_from = "year",
values_from = "pop"
)
# pivot longer to return to format of solution_original
solution_long <- solution_wide %>%
pivot_longer(`2010`:`2015`, names_to = "year", values_to = "pop")Joining Tables
Joining Tables
Oftentimes we want to use data that is split into two or more data.frames. This section will go over the use of functions to “join” or combine data from two tables into one.
We will:
Use
bind_rows()andbind_columns()for simple casesUsed dplyr commands for more complicated joins
left_join()andright join())full join()inner_join()
- Use filtering joins to filter data in one data.frame based on what is in another data.frame
semi_join()anti_join()
If you’re starting with this section, run the following code chunk to get caught up.
# call packages
library("readr")
library("janitor")
library("tidyverse")
# get data
park_visits <- read_csv("https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2019/2019-09-17/national_parks.csv")
state_pop <- read_csv("https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2019/2019-09-17/state_pop.csv")
gas_price <- read_csv("https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2019/2019-09-17/gas_price.csv")
# cut some columns from park_visits data and make sure year is nuemric data only
park_visits <- park_visits %>%
select(year, parkname, unit_code, state, visitors) %>%
filter(year != "Total") %>%
mutate(year = as.integer(year))Binding Rows and Columns
bind_rows() and bind_columns() simply merge two data.frames directly.
bind_rows()
bind_rows() adds the rows in the second data.frame to the bottom of the first data.frame. This generally only makes sense if they have many identical columns. Make sure that columns with the same data have the same name!
# we'll first make two subsets of the park_visits data
SHEN_visits <- park_visits %>%
# get data from Shenandoah.
filter(unit_code == "SHEN")
SHEN_visits## # A tibble: 81 x 5
## year parkname unit_code state visitors
## <int> <chr> <chr> <chr> <dbl>
## 1 1936 <NA> SHEN VA 694098
## 2 2016 <NA> SHEN VA 1437341
## 3 2015 <NA> SHEN VA 1321873
## 4 2014 <NA> SHEN VA 1255321
## 5 2013 <NA> SHEN VA 1136505
## 6 2012 <NA> SHEN VA 1210200
## 7 2011 <NA> SHEN VA 1209883
## 8 2010 <NA> SHEN VA 1253386
## 9 2009 <NA> SHEN VA 1120981
## 10 2008 <NA> SHEN VA 1075878
## # ... with 71 more rowsACAD_visits <- park_visits %>%
# get data from Acadia.
filter(unit_code == "ACAD")
ACAD_visits## # A tibble: 98 x 5
## year parkname unit_code state visitors
## <int> <chr> <chr> <chr> <dbl>
## 1 1919 <NA> ACAD ME 64000
## 2 2016 <NA> ACAD ME 3303393
## 3 2015 <NA> ACAD ME 2811184
## 4 2014 <NA> ACAD ME 2563129
## 5 2013 <NA> ACAD ME 2254922
## 6 2012 <NA> ACAD ME 2431052
## 7 2011 <NA> ACAD ME 2374645
## 8 2010 <NA> ACAD ME 2504208
## 9 2009 <NA> ACAD ME 2227698
## 10 2008 <NA> ACAD ME 2075857
## # ... with 88 more rows# now we will combine the 2 data.frames into one.
SHEN_ACAD_visits <- bind_rows(SHEN_visits, ACAD_visits)
# check that data from both parks is there
SHEN_ACAD_visits %>%
pull(unit_code) %>%
table()## .
## ACAD SHEN
## 98 81rm(SHEN_ACAD_visits, SHEN_visits, ACAD_visits)You generally would not split up a data.frame and then put it together again like this, but it can be useful if you have data that has the same data structure but starts out broken up by site, date, park etc. and you wish to combine it into one data.frame.
bind_columns()
bind_columns() adds the columns from the second data.frame to the right of the first data.frame.
# we'll first make two subsets of the park_visits data
visits_left <- park_visits %>%
# get first 2 columns.
select(year, parkname)
visits_left## # A tibble: 21,174 x 2
## year parkname
## <int> <chr>
## 1 1904 Crater Lake
## 2 1941 Lake Roosevelt
## 3 1961 Lewis and Clark
## 4 1935 Olympic
## 5 1982 Santa Monica Mountains
## 6 1919 <NA>
## 7 1969 <NA>
## 8 1967 <NA>
## 9 1944 <NA>
## 10 1989 <NA>
## # ... with 21,164 more rowsvisits_right <- park_visits %>%
# get last 3 columns.
select(unit_code, state, visitors)
visits_right## # A tibble: 21,174 x 3
## unit_code state visitors
## <chr> <chr> <dbl>
## 1 CRLA OR 1500
## 2 LARO WA 0
## 3 LEWI WA 69000
## 4 OLYM WA 2200
## 5 SAMO CA 468144
## 6 ACAD ME 64000
## 7 AMIS TX 448000
## 8 ASIS MD 738700
## 9 BIBE TX 1409
## 10 BICY FL 81157
## # ... with 21,164 more rowsvisits_bound <- bind_cols(visits_left, visits_right)
visits_bound## # A tibble: 21,174 x 5
## year parkname unit_code state visitors
## <int> <chr> <chr> <chr> <dbl>
## 1 1904 Crater Lake CRLA OR 1500
## 2 1941 Lake Roosevelt LARO WA 0
## 3 1961 Lewis and Clark LEWI WA 69000
## 4 1935 Olympic OLYM WA 2200
## 5 1982 Santa Monica Mountains SAMO CA 468144
## 6 1919 <NA> ACAD ME 64000
## 7 1969 <NA> AMIS TX 448000
## 8 1967 <NA> ASIS MD 738700
## 9 1944 <NA> BIBE TX 1409
## 10 1989 <NA> BICY FL 81157
## # ... with 21,164 more rowsrm(visits_left, visits_right, visits_bound)bind_columns() assumes that the rows are in the same order in each data.frame, and that there are no rows in one data.frame that are missing in the other. Unless you are sure that these are true, you should avoid using this function.
bind_rows() challenge
Using the park_visits data, make 2 smaller data.frames:
- Atlantic will have just data from Puerto Rico (PR) and the Virign Islands (VI).
- Pacific will just have data from American Samoa (AS), Guam (GU) and Hawaii (HI)
Then, using bind_rows() combine the 2 data.frames.
bind_rows Solution
# Atlantic islands
Atlantic <- park_visits %>%
filter(state %in% c("PR", "VI"))
# Pacific Islands
Pacific <- park_visits %>%
filter(state %in% c("AS", "GU", "HI"))
# All islands
Islands <- bind_rows(Atlantic, Pacific)
rm(Atlantic, Pacific, Islands)
Basic Joins - left_join() and right_join()
left_join() and right_join() are a safer way to join columns from one data.frame to the columns in another.
For either of these functions, you need to specify one or more key columns using the by argument. Key columns are the data that is used to determine which rows in once data.frame match the rows in the other. The order of the rows in the data.frames does not matter.
The syntax is left_join(x,y, by="name") where x and y are data.frames.
The only difference between left_join() and right_join() is what happens when the rows in the data.frames don’t match up one to one. left_join() keeps all the rows in x and drops any rows in y that don’t match, while right_join() keeps all the rows in y and drop any rows in x that don’t match.
Let’s combine the gas_price data with the state_pop data. The gas price data does not vary by state, so we just have to match the data by year. In order to make the way these functions work a little clearer, we will first filter the state_pop data to only include data before 2000, so that there is less of an overlap between the 2 data sources.
# filter the state data
state_pop2 <- state_pop %>% filter(year < 2000)
# left_join example
state_pop_gas1 <- left_join(
x = state_pop2, # first data.frame
y = gas_price, # second data.frame
by = "year" # key column
)
state_pop_gas1## # A tibble: 5,100 x 5
## year state pop gas_current gas_constant
## <dbl> <chr> <dbl> <dbl> <dbl>
## 1 1900 AL 1830000 NA NA
## 2 1901 AL 1907000 NA NA
## 3 1902 AL 1935000 NA NA
## 4 1903 AL 1957000 NA NA
## 5 1904 AL 1978000 NA NA
## 6 1905 AL 2012000 NA NA
## 7 1906 AL 2045000 NA NA
## 8 1907 AL 2058000 NA NA
## 9 1908 AL 2070000 NA NA
## 10 1909 AL 2108000 NA NA
## # ... with 5,090 more rowsA couple of thing to notice. The state_pop_gas1 data.frame has 5100 rows because that is what is in the state_pop2 data.frame. Every year is repeated 51 times (once for each state + DC) so the gas price data for a given year gets repeated 51 times as well. The gas_price data only goes back to 1929, so for years before that, the data just has NA for gas prices.
Here is the exact same thing, but this time using right_join()
# left_join example
state_pop_gas2 <- right_join(
x = state_pop2, # first data.frame
y = gas_price, # second data.frame
by = "year" # key column
)
state_pop_gas2## # A tibble: 3,637 x 5
## year state pop gas_current gas_constant
## <dbl> <chr> <dbl> <dbl> <dbl>
## 1 1929 AL 2644000 0.21 2.38
## 2 1930 AL 2647000 0.2 2.3
## 3 1931 AL 2649000 0.17 2.18
## 4 1932 AL 2653000 0.18 2.61
## 5 1933 AL 2661000 0.18 2.66
## 6 1934 AL 2685000 0.19 2.67
## 7 1935 AL 2719000 0.19 2.62
## 8 1936 AL 2743000 0.19 2.67
## 9 1937 AL 2762000 0.2 2.63
## 10 1938 AL 2787000 0.2 2.64
## # ... with 3,627 more rowsNow the data is very different. It has 3637 rows, which seems pretty random. What has happened is:
- All data pre-1929 was dropped.
- Data for the 71 years from 1929 to 1999 has been repeated 51 times to match the 51 states (=3621 records).
- The years 2000 to 2015 (15 records) each appear once as there is no matching data in the
state_pop2data.frame.
So, all the data from gas_price is present, if necessary it was repeated to match up with the state_pop2 data.
Basic Joins - inner_join() and full_join()
In the previous section we used left and right joins to get all of the data from one data.frame with matching data from a second data.frame. Now we will look at 2 other joins:
inner_join()- only rows that have matches in bothdata.framesare keptfull_join()- you get all the rows from both datasets even if some don’t match.
inner join:
# inner_join example
state_pop_gas3 <- inner_join(
x = state_pop2, # first data.frame
y = gas_price, # second data.frame
by = "year" # key column
)
state_pop_gas3## # A tibble: 3,621 x 5
## year state pop gas_current gas_constant
## <dbl> <chr> <dbl> <dbl> <dbl>
## 1 1929 AL 2644000 0.21 2.38
## 2 1930 AL 2647000 0.2 2.3
## 3 1931 AL 2649000 0.17 2.18
## 4 1932 AL 2653000 0.18 2.61
## 5 1933 AL 2661000 0.18 2.66
## 6 1934 AL 2685000 0.19 2.67
## 7 1935 AL 2719000 0.19 2.62
## 8 1936 AL 2743000 0.19 2.67
## 9 1937 AL 2762000 0.2 2.63
## 10 1938 AL 2787000 0.2 2.64
## # ... with 3,611 more rowsNow the data.frame has 3621 rows = there are 71 years of data from 51 states. This is the smallest data.frame from these joins. The years in state_pop_gas3 are only those that are found in both data.frames.
full join:
# full_join example
state_pop_gas4 <- full_join(
x = state_pop2, # first data.frame
y = gas_price, # second data.frame
by = "year" # key column
)
state_pop_gas4## # A tibble: 5,116 x 5
## year state pop gas_current gas_constant
## <dbl> <chr> <dbl> <dbl> <dbl>
## 1 1900 AL 1830000 NA NA
## 2 1901 AL 1907000 NA NA
## 3 1902 AL 1935000 NA NA
## 4 1903 AL 1957000 NA NA
## 5 1904 AL 1978000 NA NA
## 6 1905 AL 2012000 NA NA
## 7 1906 AL 2045000 NA NA
## 8 1907 AL 2058000 NA NA
## 9 1908 AL 2070000 NA NA
## 10 1909 AL 2108000 NA NA
## # ... with 5,106 more rowsrm(state_pop_gas1, state_pop_gas2, state_pop_gas3, state_pop_gas4)This data.frame has 5116 rows, there most of all the data.frames. This includes all data from either data.frame regardless of if there is a match.
Basic Joins Challenge
Using the park_visits and the state_pop data create a new data.frame with all the data from park_visits and any matching data from state_pop.
Hint if you need to use 2 columns as keys the you write the by= argument like so: by=c("Column1","Column2").
Basic Join Solution
park_visits_pop <- left_join(park_visits, state_pop, by = c("state", "year"))
rm(park_visits_pop)Filtering Joins
Sometimes you don’t want to combine two tables, but you do want to filter one table so that it only has data that can match data in another table. semi_join() is used for this
semi join:
# semi_join example
state_pop_gas5 <- semi_join(
x = state_pop2, # first data.frame
y = gas_price, # second data.frame
by = "year" # key column
)
state_pop_gas5## # A tibble: 3,621 x 3
## year state pop
## <dbl> <chr> <dbl>
## 1 1929 AL 2644000
## 2 1930 AL 2647000
## 3 1931 AL 2649000
## 4 1932 AL 2653000
## 5 1933 AL 2661000
## 6 1934 AL 2685000
## 7 1935 AL 2719000
## 8 1936 AL 2743000
## 9 1937 AL 2762000
## 10 1938 AL 2787000
## # ... with 3,611 more rowsThe columns in state_pop_gas4 are the same that are in state_pop2. However,all the data from year prior to 1929 has now been removed, because pre-1929 data is not in the gas_price data.
anti_join() is the opposite of semi_join(). It will keep the data one data.frame that does not have a match in the second data.frame.
anti_join:
# anti_join example
state_pop_gas6 <- anti_join(
x = state_pop2, # first data.frame
y = gas_price, # second data.frame
by = "year" # key column
)
state_pop_gas6## # A tibble: 1,479 x 3
## year state pop
## <dbl> <chr> <dbl>
## 1 1900 AL 1830000
## 2 1901 AL 1907000
## 3 1902 AL 1935000
## 4 1903 AL 1957000
## 5 1904 AL 1978000
## 6 1905 AL 2012000
## 7 1906 AL 2045000
## 8 1907 AL 2058000
## 9 1908 AL 2070000
## 10 1909 AL 2108000
## # ... with 1,469 more rows# clean up
rm(state_pop2, state_pop_gas1, state_pop_gas2, state_pop_gas3, state_pop_gas4, state_pop_gas5, state_pop_gas6)state_pop_gas_6 only has 1479 rows, because it only contains data from before 1929 when the gas_price data starts.
Filtering Joins Challenge
Filter the park_visits data so that it only contains data that can be joined with the data in and the state_pop data create a new data.frame with all the data from park_visit and any matching data from state_pop and gas_price.
Note that the state_pop data has some NAs for population for Alaska and Hawaii prior to 1949. For a bigger challenge, don’t include data those state / year combinations in your final data set.
Hint: the function is.na(data) will be TRUE when a value is NA. If you put an exclamation mark in front !is.na(data) it will tell you when the values are not NA.
Filtering Joins Solution
## In this case we filter out the NA values in the semi_join function. You could also filter them our first and save that to a new data.frame and then do the join.
park_visits_filtered<-semi_join(park_visits, state_pop %>% filter(!is.na(pop)), by=c("state", "year")) %>%
semi_join(gas_price, by="year")
park_visits_filtered## # A tibble: 19,964 x 5
## year parkname unit_code state visitors
## <int> <chr> <chr> <chr> <dbl>
## 1 1941 Lake Roosevelt LARO WA 0
## 2 1961 Lewis and Clark LEWI WA 69000
## 3 1935 Olympic OLYM WA 2200
## 4 1982 Santa Monica Mountains SAMO CA 468144
## 5 1969 <NA> AMIS TX 448000
## 6 1967 <NA> ASIS MD 738700
## 7 1944 <NA> BIBE TX 1409
## 8 1989 <NA> BICY FL 81157
## 9 1988 <NA> BISO TN 613011
## 10 1941 <NA> BLRI NC 895874
## # ... with 19,954 more rowsrm(park_visits_filtered)Day 3: Data Viz
Intro to Data Visualizations
Goals of this Day
- Learn rules for effective data visualizations
- Create simple plots, using base R functions
- Create customized and interactive plots, using ggplot functions
- Arrange multiple plots on a page, using cowplot and patchwork functions
- Create interactive maps, using sf and tmap functions
The Power of (Good) Data Visualizations
Data are useful only when used.
Data are used only when understood.
Consider three examples…
Example 1
Most people can understand this…
![Average daily Covid cases per 100k people, by region (Sources: State and local health agencies [cases]; Census Bureau [population data]](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAABBwAAAHPCAIAAAB/a150AAAgAElEQVR4nOzdeVxU5f4H8O+ZhQFmhmVg2ERcUERBRUFNDdHEMgXLq5ZaWZmit252u2333iyX6pal3l92uxlmtzSlRTMFbVFLIXEDRUMFFBVEtoFhmYFh9t8fB4ZhG5YZGNDP+9Uf8zzPec55zhl8db7zbIzRaCQAAAAAAICu4ti7AQAAAAAA0LchqAAAAAAAAKsgqAAAAAAAAKsgqAAAAAAAAKsgqAAAAAAAAKsgqAAAAAAAAKsgqAAAAAAAAKsgqAAAAAAAAKsgqAAAAAAAAKsgqAAAAAAAAKsgqAAAAAAAAKsgqAAAAAAAAKsgqAAAAAAAAKsgqAAAAAAAAKsgqAAAAAAAAKsgqAAAAAAAAKsgqAAAAAAAAKsgqAAAAAAAAKsgqAAAAAAAAKvwbHgu5WeftHuMaNmfm+UUH1gZsz5t+WdpK8IaM9O2BK786oWvzqwKNmUZM94f/9SRJbt+eT7ERu3tqJ9vvNuRwx4Y9I+WmWlbAlfucHeavPbAhzHuTYsMpXufn/3OaSMnes2R92LdtBkfTFxe+PbBf8/0anmasu9WRW0b+m3P33vH7Tl/uyOHzR/Tr/UCY8b745/61lj/Byl58kub32zxgZUxbzm18YQ7LdWXafeYSUXGVnIb7rTZ37zFP4DOse2dAgAAALTL/j0VXvfMmMAYbtyqbMwyZvy6U0zGGzdLG/O0Fw5/S9P+9mgXXjTLvlsV8vfEyvYP7DaqE2u/zWieWXrq8Gmj/Z9/A3s+JW3GBxHjlsnfPJaWlpaWlnbm0D+G7ngiYtyLP5W2X7ed047fkmWjRtrctjibtc1+d2r/f1wAAADQG9j/pZbjFTV1ovrEkd8rGnIMstw84pAx5djpxpeV8vxcYgYN7IO/vEqWPPEIo9ubcqlpdlnKkZSoJY9PYAxsmh/2StrZu/On5bIfPv/CafLa12Ld2DTHa95HB1+fQMc3f3PJcs2+6/Elj5Nhxxd4HQcAAIA7gv2DCiLPQUPqVCevlzSkS08dvjjxqUWTNOaRRsHNw5Inpge3fgaLjAU3TvJt0tAuYqLve0Ih33nU/IdkQ+nxY6dm3DdloN1a1Yz9npKh9PixkwKnIYPMh4dxvOZ9fPaclSOgyvNzrWxb93Ga8tetSyqOrF9nZW8My253avd/XAAAANA79IaggkZPWUmGHccbBggV3DzsNHT6Y9GTGiMNY8avO8VjB9YPx0/bEhjR4FOzYUXajA9M+fd/dIkaxtV8a+QdWRdtyux5zW6Q6gOn6Ht8DaYcbcYH5gN+DKV7nxs3lr2Xvx8objzG/JaNGe+PCzMbfNJkLEpfeUocaeAAMjQLupozZrw/Lqy+0WZPqfjAyqbDfsq+WxVy/0eX2A8x69PIsOPxplXI7MlEvpBUQXYzduH7lnpj2rhl098J+z3+/cD1Ltxpy7+BBmXfrQppvci8PeO3ZNn7zwYAAAB6lV4RVHD9BjdOqzBm/LpTPC8yxOueGROMX7Av4gZZbh4zdeoEN/al58XcDw83DL6/uHws+8ZcfGDlxOWFbx9KS0tLSzv7WfSOxyJfSFKGvZJ29rNHGF30miNpaWn2mujMHz2j6QiospQjKZOj73VnWpvI2zCHW/zmr+wcg7/SfzakOhERP+zp1yapTOfRXjj8rZF37ubthlrHj52a0feeEhO25I176t+JW5sYYCjd+9z4paYZF0lvqFbPahImtcZzwZZLSW9GEGfJV2lpjePKjCmrZ0WkT8llbz82dfWj9nsV5njNe+ON8fIdn7XsrGj3ln98O2zism+IiBhJZ++0yd9A2k/Lrz7SEA+Ufbcqaguzgf2zYf882KJm7dn6+IePj3vxqF+v+McFAAAAvUGvCCrMp1Ww8cMAv/ofsNlIo/TU4dM0eKAXFR9YveHk/a+/Xr+SEvta9tX/kirY/o1J0yewb1RM2KtnM46/PtGON9UEE2Y+Asrs7b91575+9TTz1FMNcwx85rz92iQVEbFDxUznuZC8dfKkyaYk2/sxoQ8+JZ85W09+9igR1YcWTXtXzn396sVJ619r+jTY2+k0JvLtQw1rLjFh9z2hUF27YcfOCp/YlY/Qb++80/xe2rllY8qJU0u/SktLS0t7L7aNv6I27tRQuvett85Ev7mmYfaO57zV64fu3PZTaf04tMefblimzKyWvvD6aYoy/cVGrMpNO/v6OBs+CAAAAOjjekVQYT6twvRmzL7TsJFGw4SKspQjKY3vxERE5BEQqEo9erqU/AfOUJ1YO8dsmAfHy8O9tYvZhfkIqMZ7bJUx49ed4rYmkJidp+xG7qwHX3/mEePnbLK+96NvPiV+2CumpZ/YyevblkV8mlH/NJrOuPCMjI5kb6fPY8Je3PZY88XBOnDLy+NXdWV+UX18HmUe0HKkgQOMx46drmTnsawIo+IDK9m4buWO+iZw/QZPoOOrZ5uPrfL0uBsXFQAAAIDW9ZKggkZPWUnGGzdL5Q1vxkRE/gNnqE5eLzFeMJ9QoTqxdkZEo/pBIEQ+c7YmvRnRWGr1gqS2ZTYCqizlSErj78FdOs+NW5WG0uPHTg0d6D3qvicUe1MukbHgxqlo08tiH31K1DBFmw0tutgd0aewQ9q2xW25Ymx/4wvbMKasnmX2xzFumWl7EDaciHnLiR0ctXVJ/eM3rcdlqohlZAEAAMCcLTe/swbXb/AESjh+apr7Sf6gpfVvxl73zJhgfOdI4sA8ZmrsBDeiMrK4LZrPnK1pc4ga5iSsnv0iHfz3TGlP3YNlTNiSN+759u2jF+/VbDi14qstXT/PfU8oXjzy+wXj0auPPxdMVDxwhurojdILVxInzlje8ONx33pKpRkX+WGjmq7+FDV14psX7dainuT58NKnNqTu2P6/yT0U5HOWNNlWsoE244OW21A2VvKa9/HZeUTETr3YsC7673TkvZhubioAAAD0Eb2lp4KdVnF4/fpvmaVRYU0yv1i/vmGwkGdkdGTzZYKMGfsTK4jIUJpxqdRUcd7/bVtExpS8wp68iXawU8+XLv+6nbVx2xvrP3rKStWJtcvXn5kXGcKedtTJI+98vquhh6ePPaW0LYGzli3d0OyXb2PBjZN8pyGD3Ft5Gq2M7+rT+GGvbF1ScSL1RH26O2/ZfP0Dk/MHfqhg16XlLIlqLaIgKruYIW/47Lngw12PMDrTCgEAAAAAvSWoYKdVEJHTxMHeLTMbxpf7xK58xPj58sYpAWXfvfDYW+tm/Pe3/c/PXvpkTOPaQReSt9a/ITH+gyZqe8MLEBsjEREbDFjAhg2m9+ziA6vZ1Z9Y7Ago0/sfxytq6j2HT5ycZhr71LeeUsSq069NUh1ZF91kbdwXHvuWWbrl+RAiGrvw/VGpbzZ7Guz4MY+AQLO1esu+WxVl/qA8AgKb7cvea41d+L5pG0SyeMut6vidcrzmPf1E1bblTZblXb7+7Uc/utT0YZI244OGORXl362KWrrsftPUeXbZsXmRIb3nHxcAAADYV28Z/kTstIod35gmVJhnNr6FM2Gvnj0yaFXUjIi1bEb0mhNpW9yIiM4+lLYl8PGIHWy+0+QPD59hX8I8H1761IZlT0Z8aWlQUI/wjIyO3HBqRBs/Bjfih71y8jOauCw6Yh0RUfSbX7w26dFtpmIm7L4nFIm5pujLMzI6cgtj9jN2H3tKngu2XHo44wPT/RKR5Mlv07bUt4HjNe/jM4Hvj58asY5HVL+uEbt4ET/slaQ3c2OWRWxj8w+mvPZ2pOlB8cOefm1i1OpZEauZyLcP/ru9p25PHK95b7xxOOatxmRbt9yqTt1pxKrcpIErY2ZFrCYiIqfJaw+nbXUnIgppfJhEkie/THozN+ato6dLYxq+oPois/b0nn9cAAAAYE+M0dj6VgkAAAAAAAAd0XuGPwEAAAAAQJ+EoAIAAAAAAKyCoAIAAAAAAKyCoAIAAAAAAKyCoAIAAAAAAKyCoAIAAAAAAKyCoAIAAAAAAKyCoAIAAAAAAKyCoAIAAAAAAKyCoAIAAAAAAKyCoAIAAAAAAKyCoAIAAAAAAKyCoAIAAAAAAKyCoAIAAAAAAKyCoAIAAAAAAKyCoAIAAAAAAKyCoAIAAAAAAKyCoAIAAAAAAKyCoAIAAAAAAKyCoAIAAAAAAKyCoAIAAAAAAKyCoAIAAAAAAKyCoAIAAAAAAKyCoAIAAAAAAKyCoAIAAAAAAKyCoAIAAAAAAKyCoAIAAAAAAKyCoAIAAAAAAKyCoAIAAAAAAKyCoAIAAAAAAKyCoAIAAAAAAKyCoAIAAAAAAKyCoAIAAAAAAKyCoAIAAAAAAKyCoAIAAAAAAKyCoAIAAAAAAKyCoAIAAAAAAKzC675TG2QpX++7UiMOnbdokqQx+/L++JQSYoiIqN/0uNmBTSpZLr3r6LL2f55cEjAlbmZwsxJZcsL3WQr2QbmEzV04Xtrx0k6dCt8XAAAAALSDMRqN3XPm+ldVxiyoYMMMfuj8BZMkRHTp0NYTBW6mt1jLpXefxvf1FkHF5f3xKWX+056ZFURExam7DmTWmB1jubSVq7R1ML4vAAAAAOiI7hr+VJz6SzaNHO5vMMuT/X7kco04dHpDv0XIrHnB4qoLR1Ll7ZfebWTJCSmCKSuWzR0ualF26VByCflPnRXEJn0mPTbZ35CffDC3A6WdORW+LwAAAADokG4Z/qTL2n8gkxs2d5Lw7EVTpkGWVaBg3EKDzYZCSYMGiLIy86/JJkWQpdKWP36z44LYz+adIeb5zQbzsD/Ds5/dG35fJyIL43/artKkGQFTotTJx9heBa7/1Gca3tGbtces6PL++BRN6IIFk4r3x6eUikc2HSFGRNIpi1YQkUHW/MaJLl8rYLj+Q80HGg0b7Hui4PbVLAoMtlzaiVMN8rDl9wUAAAAAd7DuCCouH0wudg9dMF5Kl8xyDeVyJbmEDW3y8uwpERGVyMvJQJZKqelLatOZBrLkhO/3JtC8RZNcms5AKE7ddWDf1zR34XgpOzin//S4xwLrx+3s2V7JvuLXj/+JWxRUf+b2qzSXn3xtetyKh+oPO7b9EJmNJuKGzY0bL61vp6mIiHR5P8dnKogYpjMP1yArryGm3+AmzeB4SERUopDLLZcSSTpxKovfSKe+LwAAAAC4s9l++NOlQ8ll/tNa/qhfJldaqGW5tKnLB5OLuf5TGyYJSCeN8TEq8q/JKPt6ESMOHd/we7zX0AARVcvLiUhWXtX4kzxHGrk4Lu7pcRJ5i/E/vOCHJvtX/nE2x0KVlgKmzDYddl+oUF9w7oyMDLKUXzNr3EPvb/jZXnpv9AingvQzDZ0PCoUiYEpcXFzc8ubdFJYYyi09KculnTrYdt8XAAAAANzhbNxTUZy660RB/+lxrfycb4FCbmkcfrOf2HVZV0uICTD7fZ0X/FBcsKxMRp6zVoYQNV2SiK0p9XA1ZhUcM+8o4Eg9JXQ5pcX4Hw83sT7zai7NbqNKO0y/1pfK85v9ls+Regjpct5VORtmmMVFtlFdWUZtt6+6sowsFLc82K3N0k59XwAAAABwx7NlUMH+Nh8w5bFW1xX1lIiIVK1WFEsknmSptAMXl3pKyRROcP2nxc0KYpcnYotDZq2kQ1tPFByLjz9GRObrn+obM03cLFexwGysERFVZ+yLz2hSzrh34GbaPX+rRS5unpZLO3Wqbv6+AAAAAODOYcug4srZS0riKJPj45NNeRyizD3xmYw4dG60RERX8q7Kx0sbXzrL5EoiF4kHcchSacfIkhNSSsUj57cxmiikvh+DnZJx+2j8QVo+mNqYgd1mlfbiCnZMUYBEQnJqOw4p7uAttcSRegjJePt6DgU39tWYLmq5tHOn8uju7wsAAADgbpSTkyMSifz8/OzdEFuy5ZyKkFkr45qa7G9gxKHz4+KWL5rkKQ32Fxsr87LMhs7IcvKUjDhgiJQ4FkvN8YKHepOx6Qicy/vjP/0pS16hILcBwa3GB7qsK6Y1VXnBD80JFRIpysuGD/FvdlEyyFJOZrVdpZXlmJowvVh7DQ0Q0e2rWU1KL6eesHrF1RFD/I36gqvmS8RmXy8i6jc0uN3STpzK8jfS8e8LAAAAAEw2b968ePHixx57LCEhwd5tsaXu2qeiNdJ7o0cIFZlHG3YyuHRob5bCdXQ0201gudTciAmhoorMPT81vK+zU8NnBkvcxWT2miv7/chldjaxLmv/58kpvyaYtlBofP0dPi5EqMjcdyin4eSXE/dd/iP569SziW1VaenW+frD2AFgXP+x46XEkUaO9jfkJ39tmpldnLrr98xLZtfqopBZU7yp4FjDeYpTd50o4Jgmi1suLU7dFR8fv72h1OLBtvq+AAAAAIAUCsXLL7+8e/du9vOmTZtWrFihUCjs3S7b6L4dtYmILh3amlo1qukmDOazqFuODrJc2sjC/g8N1V3C5oYV7TvesF20+WYUzXaTaDKx22zzaQtVmjTDPzSkMvOSsrXDzHe6MNtPo8k+1q1pcukW1cnC3hqWS9n2NG2k5VPZ5vsCAAAAuJspFIoVK1bk5OQQUVBQkEKhKCoqIqKNGzdOnTrV3q2zge4NKu54TXfM6BNkyQmHHaMXY386AAAAgB6Tk5OzePFiIoqKilq7di0RJSQkZGdnr127ViwW27t1NtAtO2pDr3Xp0N6SAY8sQEQBAAAA0IOCgoI2btyoUChiY2PZnLi4OPMDCgsLjx8/HhMT00djDPRUWKUP9lQAAAAAQK8zZ86cwsJCPz+/NWvWhIeHt3pMUlLSsWPHcnJyCgsLxWLxgQMHek8EwmX7X6BrOJ7B4eHhQ5rvAAEAAAAAd7XCwsIVK1a8++674eHhHVk9Ni0tLS8vT6FQJCUlKZVKgUCQnp6uUChMdXNycp5//vmbN2+yc7s1Go2Hh8fIkSO79zY6zJY9FYWFhewH9ubZJMMwvr6+lpNGo5GdqmLDJIfD8fHxsZw0GAzFxcW2TXK5XG9vb8tJvV5fUlJi2ySPx/Py8rKc1Ol0paWltk3y+XypVGo5qdVqZTKZbZMODg6enp6WkxqNpqyszLZJgUDg4eFhOalWq8vLy22bdHR0lEgklpN1dXVyudzKpJOTk7u7e6eSKpWqoqLCyqSzs7Obm5vlZG1tbWVlpW2TQqHQ1dXVcrKmpqaqqsq2SZFI5OLiYjmpVCqrq6ttmxSLxWKx2HJSoVCw/7uyYdLFxUUkEllOVldXK5VK2yZdXV2FQqHlZFVVVU1NjW2Tbm5uzs7OlpOVlZW1tbW2Tbq7uzs5OVlOVlRUqFQqK5MSicTR0bFTSblcXldXZ2XSw8NDIBBYTpaXl6vVatsmPT09HRwcLCfLyso0Go1tk1KplM/nW07KZDKtVmvbpJeXF4/Hs5wsLS3V6XS2TXp7e3O5XMvJkpISvV5v26SPjw+Hw7GcLC4uNhgMnUpevHjxvffeUyqVWq12xYoVDz/8cLODfX19GYZplvzkk0927txZU1PD5/OJSKvVCoXCL7/8csiQIQzDFBUVPfvss0qlctSoUTk5Ofn5+d7e3l988QVbt6ioiH2rt5xk39ItJ5u90ncQhj8BAAAAANhMYmLiunXr2M+LFi166aWXOl63sLBw7dq1586dY5O+vr6JiYktD9u9e/fmzZuJ6NNPP21rrFQPQ1ABAAAAAGAb6enpK1asICKRSPTSSy+ZpmV3ChtIDBs2LCio9e0H2AnfSqVyzZo1XbuEzSGoAAAAAACwjfj4+Pj4eCLavXt3WyGBTaSnpxcWFvaSiIJ6dkdtAAAAAIDeLm1LYERERERExKcZLTLHb8lqyDGU7n1u3Nhmh8XFPfwKRyMSibo1oiCi8PDw3hNREIIKAAAAAABzYxe+P4ExENGNW5X1WcaMX3eKiYiMN26W1ufpC6+fNnKIsyRq9IX3x4VFRETc/9ElIicisvvrving+XtiZftH2wKCCgAAAACARhxp4AAyENGJI79XEBGRQZabx742G1OOna5/TS/PzyUip4mDvcl8NoH40TMXOzU520rnzp1jF9OzL+yoDQAAAABghgm77wnFtzvcVSevlxC5mzoliIjo3M3bRG5EZSlHUoicJkff6864vXo241V7tPSna+qZY8e2zOd4zfv47LyebAl6KgAAAAAAmvAfOIOocbDTheStRDR50mQiUl27UUFExoIbJ/lENKi/GxkzzIY/ERFpMz6IMBn34k+ljWcuPrDSVGI6nohMJ2k2ncN85oapbuQLSWwXyswhglav1Wz4EzshJPKFpAqzq5hPBbEeggoAAAAAgCa87pkxgTE0DHYqu3HNkZjI+6IHEJEq9ejp0oYBUZwlUWHN6xYfWDlx2TeNaWPK6ln1b/BpWwJj1qeZSuRfPsmGB4bSvc+NX/qtsXEM0bZlzadDbFsWYaqrOrH20Y8uWb5WS6oTa2eMW2a6yrblTaIdKyGoAAAAAABowjSt4tzN24bS48dOCpwmTZ8Se98jjI6MKXmFVHrq8Gkjx2niYO9mNY0ZO946RUTLP0tjJb256O1DafP9yrUZH6zc4U5E0WuOpKWlnTn0jwmMQXVi7bcZFXvffpOd8/1VWlpaWtrWJRVEdGT9OvOXfqfJaw+npaWd/ewRRkdsh0nb12rrvtgjk96MICL2Rmz2xGx2JgAAAACAOwMTdt8TCiJSXbtRVnj9tJHjNGSQO+M/aKKWiPamXCq4eZiIJkff6968Yv0x25bVj0fKuOelmV7k4eXBjqEiJnLqBDeqn/ZwLi0tbbnfr8dOCogo+o0lwURENHrKSqLmL/3112o4PxER07+ta7V+Uw39Kh4BgTZ4RM3ObfMzAgAAAAD0dey0ClXq0Xc+30VE8yJDiDwjoyOJSHV12/92uhI7oaI5z4eXPtWYMqasntV0WgUzaKBXO5fm+g1uvqZt6zzauVYPQlABAAAAANCcaVrFiVQeMZED/IgafuNXpf5ev0NFiwkVRMQPe6VhMFJEfZYxZfM3DXOyzXa6aItpsanWgpbOXKsHIagAAAAAAGjONK2CqLF7wdSHQESSJ6YHt1ZRm/EBu/G2z5ytpgkSqms3AhoGNTXsdFH23aqQiIiIf54aM3WimoiOvLWD3a7bNFCKjWQsaOtaFUbGilvvCuxTAQAAAADQQsNuFWQWP3C8oqZOfPN0qhMRjR3Yr2UlQ+nevy5PICPn8YgdDXnuRPT40zFeYTFbl2xducP9yLroiHVskZPkyS/fmzPYcM/6Y7PfOW3Y0VDLnYii31wz04sMbXdrGEq+b+ta7sxea2+/k9BTAQAAAADQivrdKprED56DhtQRNc63bobjNe/js7+8NknVmMVEvn0obUUYEVHEqtzGcUpEkie//OX5kPpaZz5nl3ViLf8s7b3YdsY+cbz/ZOFaPYwxGo3tHwUAAAAAAO2r1pHLXTgWCD0VAAAAAAA2kJ6enph43I4RRXx8/Lp16+xyafRUAAAAAAD0eTk5OYsXLyaiNWvWxMbG9vDV0VMBAAAAANDnBQUF+fr6ElFSUlLPXx1BBQAAAADAnYDtoEhPTy8sLGz3YNtCUAEAAAAAYK3ExMQVK1akp6fbsQ0xMTHsh4SEhB6+NIIKAAAAAABrJSQkpKen7969245t8PPzi4qKInuMgEJQAQAAAABgrZycHCIaNmyYfZsxdepUIlIoFImJiT153btwFV0AAAAAgDtTbGzspk2blEplenq6aQ2oc+fOtTzSx8fHz8/PVtdFUAEAAAAAYBXTVIrw8HD7toSINm3atHv3bnZ5WSLavXv35s2bWx4mFot/++03W10UQQUAAAAAwJ0jPDy8I7GNbXerQ1ABAAAAAGAVU09FUFCQfVvS0uLFi6dOnVpUVNQs37adKggqAAAAAABsQywW27sJrfDz87Ph9IlWIagAAAAAALBKeHi4r6+vaWL0XYix7WgqAAAAAAC422CfCgAAAAAAsAqCCgAAAAAAsAqCCgAAAAAAq+Tk5LA7at+1MFEbAAAAAKDrcnJy2J3mDhw40N2LLPVa6KkAAAAAAOi67Oxs9oNSqbRvS+wIQQUAAAAAQNeZ9pXrhTvf9RgEFQAAAAAAXadQKIhIJBLZuyH2hKACAAAAAKDr2OFPd3M3BSGoAAAAAACwxt08lcIEQQUAAAAAQNexi8lGRETYuyH2hKACAAAAAKDr2NkUvr6+9m6IPTFGo9HebQAAAAAA6KsKCwuLiorCw8Pt3RB7QlABAAAAAABWwfAnAAAAAACwCoIKAAAAAACwCoIKAAAAAIAuSkxMnDNnTmJior0bYmcIKgAAAAAAuigxMbGwsBBBBYIKAAAAAIAuunr1KhH5+fnZuyF2hqACAAAAAKCLFAoFIahAUAEAAAAA0DWFhYXsh7t85ztCUAEAAAAA0DVFRUXsB/RUIKgAAAAAAOgKU0+FWCy2b0vsDkEFAAAAAEBXmHoqgoKC7NsSu+PZuwEAAAAAAH1STExMdnZ2eHi4vRtif4zRaLR3GwAAAAAAoA/D8CcAAAAAALAKggoAAAAAALAKggoAAAAAgK5ISkpKSkqydyt6BUzUBgAAAADotMLCwrVr1xJRUFAQVn+yfVBhkKV8ve+Ksj7lEjZ34Xipefnl/fEpJcQQEVG/6XGzA6njpQAAAAAAvYJpPVmFQmHflvQGNg4qilN3HcisCZgStzi4IbkvXj4lbmYwUUO8wQ9dEDdJQkSXDm09Gv91eUPUYbkUAAAAAKD3SE9PZz+gm4JsO6fCIEv5NbPGPXQ+G0IQkc+kMd5kvHU+VU5EJPv9yOUacej0SRK2NGTWvGBx1YUjHSkFAAAAAOiNsJ022bqnYvjMuEhJi1yOq6eEyCDLKlAwbqHBZgdIgwaIsjLzr8kmRZCl0padFbqs/Z8nl7CfGXHovEX1sYh5frPBV2wvCvvZPXT+gkmmS8mSE77PUjCdqQIAAAAAd7Xs7Gwi8vX1tXdDegVbBhUcqWezl+5Lh5dTlxUAACAASURBVJJLqP/0WUFEZCiXK8klbGiTQzwlIqISeTkZyFIpNQ0q2MghoH5UlSw54fu9CTRv0SSXJvns4Kuvae7C8VK6dGjriYL+0+MeC6wfZ7Vne+XUZ2YFsbM4yvynxS0Kqj9z+1UAAAAA4G7HTqVAUMHqjtWfGidbM+JR8+Pqf94vkyuJuG3VsVza7PwHk4u5/tMaBllJJ43xyUrOvyabJLxexIhHjW8YfOU1NECUeUVeTiSVlVcxXP+h7LRvjjRycVykQVYmJyoyC3uIiBf80OTrW0+dzRk/y72tKuitAAAAAIDi4mIi8vPzs3dDeoXuCCpGPBQ3gv1UnLprT3xmuwOHFHJL8yYU8iZv8rqsqyXEBAxu7DHgBT8UFywrk5HnrJUhRE2XkGJrSj1cjVkFx7YfIlNXA0fqKaHLKQWNkQPLw02sz7yaS7PbqAIAAAAAQCKRiIjCw8Pt3ZBeoXv3qfCZ9Njkyq0nMk/mTpo9QCIiUrV6mFgi8SRLpR24lNRTSqZwgus/LW5WELucFFscMmslHdp6ouBYfPwxIjJfr1bfmGniZrkKAAAAANzldu/enZOTg6WfWN2++Z2Hm5gKFOUyGuQhEdGVvKvy8dLGIKFMriRykXgQhyyVdowsOSGlVDxy/qLWu0VC6vsx2CkZt4/GH6Tlg8niDOxWqiCuAAAAAAAiwmKyZmy4pKwsOeHT+PiDuU0zc/KURGIPKXGkwf5iY2VelrxpKSMOGNJeqTle8FBvMjYdMXV5f/ynP2XJKxTkNiC41fhAl3XF1DBe8ENzQoVEivKy4UP8m12UDLKUk1ltV5F1/IEAAAAAANwVbBhUSKcsivSmgqNmcUVx6i9ZCiZgCvvrvvTe6BFCRebRhp0nLh3am6VwHR3NdhNYLjU3YkKoqCJzz09Z1HBkcpn/tJnBEncxmUUIst+PXGY39tZl7f88OeXXBNOWF43hyvBxIUJF5r5DOQ0nv5y47/IfyV+nnk1sqwoAAAAAAJhjjEajbc946dDWEwWmWKXlPATzWdSdLW1kvh8F19+00qt5dZewuWFF+46X+U97ZlZQ080ozKs0q0WmFWktVgEAAACAu1diYuK6desWLVr00ksv2bstvYLtgwoAAAAAgDtbXFzcuXPnxo4dGx8fb++29ArdPlEbAAAAAKD46935a5/TVVVaPozvFeD3/D/7LVvRM60CW7HhnAoAAAAAgFZcX7v6+ouPtRtREJG2ND/vjZVnRw+4/dmnPdCwLjt37hwRRURE2LshvQV6KgAAAACgu+hrarKWP1n1214ichwyRhI738LBRp2u7OttWlkBG1oUfvQv9Fr0FbacU1FYWMh+YLcrZ5MMw/j6+lpOGo3GoqIi2yY5HI6Pj4/lpMFgYPdXt2GSy+V6e3tbTur1+pKSEtsmeTyel5eX5aROpystLbVtks/nS6VSy0mtViuTyWybdHBw8PT0tJzUaDRlZWW2TQoEAg8PD8tJtVpdXl5u26Sjo6NEIrGcrKurk8vlViadnJzc3d07lVSpVBUVFVYmnZ2d3dzcLCdra2srKyttmxQKha6urpaTNTU1VVVVtk2KRCIXFxfLSaVSWV1dbdukWCwWi8WWkwqFQqFQ2Dbp4uIiEoksJ6urq5VKpW2Trq6uQqHQcrKqqqqmpsa2STc3N2dnZ8vJysrK2tpa2ybd3d2dnJwsJysqKlQqlZVJiUTi6OjYqaRcLq+rq7My6eHhIRAILCfLy8vVarVtk56eng4ODpaTZWVlGo3GtkmpVMrn8y0nZTKZVqu1bdLLy4vH41lOlpaW6nS6jiT1JSUX/jRDl3+FiFyiHh6+/avymppWD/b29uZyuaWlpTq1umrfHuUX/9EWXSciIuJ49uv31zf9lzxdKpfr9XrTwSUlJdYnfXx8OByO5WRxcbHBYDBPZmRkvPbaa3w+f+PGjcHBwc1KO5L09fVlGMZysqioiH1Rt2GSfUu3nGz2St9BmKgNAAAAALZXdeZ09pKZ7JAnnxWvD17zFjFMB+sadbqi3TsL/+9tTUNowfce4PfCar/Hn2T4/O5qcYelp6evWLGCiD799NPw8HB7N6dXwJwKAAAAALCx4q93XXroHjaiGPzvXYPXvt3xiIKIGB7Pb8nTEedyB7z7GV/qT0Takry8fy5PG9P/9nb7r7bEdrQSdtQ2g54KAAAAALCl62teL47/FxHx3CTDdvzoOm68lScs3r2z4IM3NcU32aSg/7CAN96Xxs6xtqFdpVAo4uPjg4KCYmNj7dWG3gZBBQAAAADYhl6pzFr2RNXxH4jIcciYkK8TBf362erkxbt33tqwWluazyadQycNenuT64R7bHV+sAaCCgAAAACwAfXt25cWzKy7kUlErlPnBn+2kysU2vYSerW68NP/Fn60Wq+sZXNc75s/aP27zoFDbHsh6CwEFQAAAABgFX1t7e1tWws3vWTQEhH5/PmNwW+u777L6Soq8jZtKNm+wZTj+eizA19f6yCVdt9FwTIEFQAAAADQRbqKioJPPy757B19TR2bE/hhgvcjC3vg0nX5+TffXS//YTub5PDJ+8/rA55/kSsSdfelExMTi4qKFi1aJBaLu/tafQWCCgAAAADoNHVJScFH/5bt2MD2ThCRU1D44I3/tX5adqcoLl64ufYfipM/skmemyRg7X98Hl3UfVcsLCycM2cOEW3cuHHq1Kndd6G+BUvKAgAAAEAn1OXl5by0Kj3Mp2R7fUThPCoyaHvSmONpPRxREJF41OiR3x8a9sWPjoNCiUhXKb/+18VFX33ZfVdkd1ImInRTmOPZuwEAAAAA0DfUZF0p2LK5fN9nphyXyNh+z7/sHjnFjq0iIo8HZno8MLN49868NSv1ytobrzzFETh6L3i0O65l2qQCQYU5BBUAAAAA0I6qtLNF8f+VJ35hynG7f2H/l/4uHjXajq1qxmfxE87DQy7PizCojLmrFnIcnbpjL4ucnBz2A3a+M4egAgAAAABaV3XqZFnS/vLvt+sqykyZHnOX9f/bq85DhtqxYW1xGTN2+O4Tl+ZOIqKrcQ9xdv7iET3Dtpcw9VSAOQQVAAAAANBEdXqabN8eeWKCaac5IuKKRe4xTwX89RXHgAA7tq1drvdMDE74NWvRfUR0den9nK+OuU+JsuH5s7OziWjs2LE2POcdAEEFAAAAABARKf+4KPthb/n+BM3tq6ZMrljkNvNx6dwF7lOiGC7Xjs3rOMnUacN2/Jzz9AMGLWUvmTp8zxnXiHH2btQdDkEFAAAAwF1NXVRU9NUX5Xu/UuddNmVyHMk9dqVnzMMe9z9gx7Z1mceM+wP/+8O1FQ8b1JS1cPKIvWfEo8NscuZhw4adO3du8eLFNjnbHQP7VAAAAADcpUr27pF9/WX170mmHI4Duc16xmPWnO6Y4tzzihO+uv63J4iIKxaF7j8tHD7C3i26YyGoAAAAALi71GRnFe/4X9l3/9UrlKZM8cQHvRY96TnnYa5AYMe22dzt/32W98/lRMRzk4QmnXYOHGLvFt2ZEFQAAAAA3BX0SmXxnm9kuz6vzUw1ZXJdXD0f+bPvU8/cwW/btz7ecuvtF4iI7+k76qezGheX48ePFxYWNjvMz88vJibGlMzJyTl27FjLs8XExPj5+XVrg/siBBUAAAAAd7jq8+eKtn1SkfQZuwE2y3lUpM/TK6Rz599hXROturnhncL/W01EDv5Bt1atfufDD1s9bPfu3abdJ+bMmdMy8CCiqKioTZs2dV9T+yhM1AYAAAC4Y5UdOnj7Pxtrzjf+4s4RkORPf/FbGicKHWnHhvWwga+9rlcoSrZv0BTkGN9f4+03tKS8vNkxvr6+vr6+puTYsWNbDSrCw8O7t619E3oqAAAAAO5Atz/fVvzpZnV+linHcfBI72WrfBYs5IpEdmyYHV19+a+yXR8SkdPwCaMOHLlrn0N3QFABAAAAcOfQyuW3P9ta+r9Nukp5fRZDrtPm+8Y9J4maatem2dnmzZtFQmHkyeMVPycQkTD8vpHf/8hxcLB3u+4QGP4EAAAAcCdQ3bx566PN8j0fGTT1OVwh3/PRF/ut/Itj//52bZr9JSYm7t69m4iiduxwqVFU/55Uk/5r1tLHR+z4mjgce7fuToCeCgAAAIC+TV1YmPvPlyt/TjDlCPoP835mle8TT3Gdne3YsN6DnXXt6+ubmJioV6sz582qSf+ViDwX/SVo80cWKhrr6oxqtVGjNqrVpK4zqNWk1nAHDORKJD3V9r4BPRUAAAAAfVjN5UuXH52hLStik65T5/r9+QX3KVH2bVWvkpiYyE65jouLIyKuQBD6beLFmPtUV06XJfyH6+Qc+M4Go0ajLynSF5cYykqNdXWkURvUatKoWz9j+mmOmzt3wCDewMFcqbQn76XXQk8FAAAAQF9VkZKc/eRUg8pIRF5L/ub/l79ipFNLbDeFSCQy33dCK5dfmDFeU5hLRL6xS6Tj7+nayRmhiDcwkDdwINf3rt68Aj0VAAAAAH1SyXffXP/bQqOOGD4N3ZrkOWu2vVvUGx07doztpli8eDER6YsK9SXF+pJifWnJ4EXPXPtkta7WUJS4g+/g6BYWRkSMUMRxl5DAkSMQkEDACAQcgYAcBIzZlG5jTY322lV9QR4RGWuU2ksXtJcukMCRFzCQN2gwL2CAne7VntBTAQAAAND35P/fpoINLxMRV+gYnHDcddx4e7eol1qxYkX62bNCHu+71a87y2Wk05qXqopLcuPXGbTEMBT43ufS+Y8yHZ6FYtRotDeu63Ov6gsLmhTw+Nz+Ax3GhnPd3W11F70fggoAAACAPsVgyP7rc+XfbSUivldAyN4jzkOG2rtNvZG+rEyXnzfxkUeNOt3TkfcujZrcWCZw5Hp5c719ud7eyhs3Li+816glhk8j9pxyHT+hsxcyqtXa3Kv6G7n6oiab5XEDgwQR4zlisfX30vshqAAAAADoM/RqddaTC6uO/0BEjoFhId//LPDysnejehOdTnf7lj4/X5d/06iqJaKMvPxzN28tjZpMPD7X14/brz+3n3+zPoSyHw/lPDObjMRx5o48+IcweHjXLm6srdXmXtPfyjPvu+APG+EwNoIRCq25rY42oK5Ol5+ny8vT386v75PhOzBcLnE4xOEwHA7Vf+YyPC5H7CqYYrOtSxBUAAAAAPQNuoqKzEdiajNTiUg0/v6QXXuxJ7Q5vUym+jGx2ZJNHKk3t19/nn9/ro+PhbpFu3bcePlJIuJJvEYePO00cKA1LTEqFJrz6dqcK6YcfshohzFjGUdHa07b5uVqarQ5WbpbeYbSkk5VdF6wmOPqapM2IKgAAAAA6ANqb9y48uiD6lvZRCR5aOmwj7YyfL69G9WLGOvqavftMdYoiF2Rqf8Arr8/18+f6fCe2Xmb3r+98TUi4nv6+v/9Pd/HlljbJIVCk3FOm325Ps3j80NHOYwK63iTLNMV3NLfLtDfvmWQlzcr4nh5M1yu5eocF3dB5BSbtIQQVAAAAAD0ftXnz11ZeJ++uoqIfJ9bM2j1Wnu3qNdR/Zikv32LiART7uMHDVMoFMePH4+KihJ3ZkrDtVf/Vrrz3+xn4egpgRv/IwodaWXDmvdaCBwdRo9xGBFKPEursBq1WtKoDXVq0miM6jqDuo7UGqO2YZZ5bU1jrNKAcXTi9h/ACwjg9utvq7il4xBUAAAAAPRq5Yd/yV7yAPt54IbP/ZY8bd/29EKa8+c06aeJiDc02DFqGhGtXbs2KSkpKipq06ZNnTpVwdaPC95/nt36g4g8Fz43+M23eFav42RUKjXn01tGAtZjYwl+8Aiut7fNT96JZiCoAAAAAOi1ZIkHrv35IaOeOHwa+vkvHtEz7N2iXkdfVKg6uJ+IOO4S54fmEY9XWFg4Z84cIoqJiVm7ttO9Opry8htrXy/f+ykZiYi4YpH/qxv7LV1OHI6VTTUqFHXpafprWV0/hWkDDSch168f168fVyKxslU2gaACAAAAoJe6/cX2vH8sIyKOEzP8m1PYjKIlo0pV+/23RlUt8fjO8x7hiF2ooZuCiA4cOODn18WNrqszzuf+dYUq+yybdAwMC9z0ieuELm683aTNNUp9ZQUZDKQ3kMFABr3RYDDqDWTQszlGg4HhcIjdd0/gyAgE9f/1+KCmjkNQAQAAANAbmeYN8yReId/9KhwRYu8W9T5GY23SfkNJERE5zniQN2AgEVnZTdHs/EVffXnrXy/pKuVshiT2qUFvvSew60Cj3snaThwAAAAAsDGj8drfX65fich7wKgfzyKiaJU67QwbUfBHjmEjCiKKj49nP8TFxVl7AYbxfeKpsadveC35G5shT/zi/Hifku++sfbMdxwEFQAAAAC9icGQFfd06ZebiMhxUOion087BgTYu029ke5WvvbCOSLiePsIGrbBLiwsZAc+RUVFdXngUzM8F5chGzaN/vWScOw0IjJoKHfVwurz52xy8jsGggoAAACA3sKo1V56/BF50pdE5DR8wqgff8dIm1YZa2rqfjtCRMR3cJr+ADEMm3/s2DH2w+LFi217ReHwEaMP/hq45WuOIxFR1uIZ6tu3bXuJPs3S+rgAAAAA0GP0tbWXH5unOPUTEQnHTgv9Lonr7GzvRvVKBoPqyE/sztmOU+9jzJ5SbGxsenq6n59feHh4d1zZe8GjXJE455nZukr5pUdnjz5yits9m2T3OZioDQAAAGB/uqqqP+Y+oLpymohcp80b/r9dHIHA3o3qpdQnU7WXLhARP2SUYOLknm9A/r83Frz/ChG5TV8wYuc3pn6SuxmGPwEAAADYmbqk5OKD97IRhSTmyZCvvkVE0RZd3g02ouBIvQQTJiYlJc2ZMyciIiI9Pb3H2hDw4svuDz5GRJVHv7vx7ls9dt3eDEEFAAAAgD3VFRT8MWty3Y1MIvJe+mrwti+s32TtTmVUKOqO/UpExbWqnbeKpk2fvnbt2sLCQiLKzs7uyZYM+2S7c+gkIir6aE3p/h968tK9E4Y/AQAAANhTxozJtZmpRNTv5Q0DXnrV3s3pvYwqleqnpMLr1/93/PefbuaZdoLz9fWNjY21wQKynaQpK7swbZS2rIjh08iDF0QjR/VwA3oVBBUAAAAAdlOw9eP8dX8hIr+X3hv48mv2bk4vptfXJv1gkJWu2rH7gqyMEQqJaOzYsbGxsbGxsfZqlPKPi5kxow0a4nv6jjpy/m5eqgtBBQAAAIB91OXnZ0QOMGjIafiEsMMnGC7X3i3qvVSHf9Hn5RLRzwVF7yYdiomJiY2N7aYlnprRGYw1ap1Cra/R6GrUOqVabzB/f0456vD6QiLSDxpriN9vFDgSkQOP4+LId3HkiQVcF0c+l3Pnz+RGUAEAAABgHxfn3K88e5jh0qgjl4XBw+3dnN6r7Lejjrk5RMSReDjP+RPxurgpQmWttkihNhqNHIbhMgyHIQ6H4TAMhyGGYTiMkcMwHA6j0hiUal2NWlej0SnUOo3OYPm0Tjs/dk94l4hqImKr1n/SygEOXLGA7+rEEwv4Lo5ckSPfkXenTZvBPhUAAAAAdlC0a4fy7GEi8nl2HSKKtigUir8tX5526tS/FvxpytgxTjNjuhBRVNdp8yvqCipqlWqdle3hchixgMfjNg0J/vw37fXL/NP7hWmJ/G9G0zMvKNW6Oq3eVK7S6FUafamisYaAzx0kcR7sKXR2uEO6p9BTAQAAANDTNDLZ+YkB+po6wYARY5LPcxrmHIO5nJycvz333O2rV4nowbBRb32xg+vm1vHqtRpdnlxVUKmqUmm7cHUBjysScIUCnkjAEznwhAKuUMB15LUeA+jr6i4+GKXKOkNEQZ8f9Hxwls5gVNRpq+v0ijptdZ1OUadVtBbSeLs4DvIQ+rv1+R30EFQAAAAA9LRLi+dX/baXGApNPOsSHmHv5vRGiYmJmzdsqLxdQEZ6cNTIVzdudAsa1pGKKq0+T64qqKitbBpLOPA4Ae7OPi6Cdmc4OHA5IgGvsxMh1CUlF6PHaMuKOA4UmtT6YlCKOm21Wqeo05cp1cXVdaZ8Jz53kKdwsKew7w6LQlABAAAA0KNK9/9wbeVcIvJ68qUh7220d3N6o82bN+/+6itDZaXRYFg1477HX3mVHxTcbq2b5bW5ZTUVtRrzTD6X09/dyd/NyUvc7fsJVp8/lzkrnIj4Ht6jjmQIfHwsHFyr0V8vr71epjSfs9HPzSlAwi+tO6w1qPgcJxcHLy7H0UXg5cRzdeK5dnf7rYGgAgAAAKDn6Kqrz40foKuq5Ev9x57M4gqF9m5Rr7Np06aE3bsNlZVCHvfdR+aNi31YMH6C5Sq3K+syi6oUdY3ji3gcpr+7cz83Rx+XHh1ZVLLn29znHyUi8cQHR35/qCNVblWobpTXlCrUphyGU8t3yHYQXGWYOvMjnXhujjwXJ56bM89VLPD2ch5q28ZbAxO1AQAAAHpO7j9f0VVVElHglh2IKFpVWFhoqK4O9JS8+8j8fqPHWI4oKlXajFuVZTX1vRMOPI6Pi6O/m7Ofa7f3S7TKe/4jFUd/kf+wXXHyx1tb/t1/1YvtVunv7tTf3alWo78mq7kpV2p0RqPBWVM3RlM3hudwk++Qw+MVs0eqdJUqXWUF5ZvqejkP8xIO9XIO4nPsc78m6KkAAAAA6CHyY79lLbqPiDz+tHzYx/H2bk4vVfbjobO/HokcNpQj9XKOeZja2L6jVqP/o7D6VkWtKWeIVDTCV+zAtfO0BH1t7fkpYZrbVxkuhR7KEI8a3fG6Z4p2lVZztJrhep3UlCkUkIe42tW5olZXpDXUKTQlOoO6WUWRg3c/0Shv4VB7jZJCUAEAAADQE/R1defGD9XKCniubmPP5PFcXOzdol4kJyfH19dXLBZrLl7QnEklIo7YxemheYxjK4OXNHrDlWJFrqzGtAmdv7vTKD/X3rM8q+LihcxZYUY9OfQbOiY5g+vs3JFa8rr8s0W7iMhPNGqg+IHcspo8ea1WXz/jgh3QFSgVujnxiai0Nqek5mppbXazAMOJ5+blPLSfeKTYoUe39+auXbu2J68HAAAAcHe68cY/qlN/IqLBHyaIR3fi1+s73rp16959993MzMwHQ0PUvx8jIhI4Os1+iCMSNTtSbzDmlNacuiGXKdVsPOHr6jhugHuQl4hv7w4KltagLq65IvD25vN9qk/8olfI6wpknrNjO1L3bNFunaGOxxGE+yx04jv4uDgOlYqEAl6tRq/WGQxGqlRpr5fVFFWrOQzjI/b2EQUNdpvoLQziMDy1voaNLnSGuip14S3F+UJlpljg3WMdF+ipAAAAAOh21elpmbHjyEiu0+aF7N5j7+b0FgqF4uWXX05PTyeisMGD/+/BGaTXEZFT7Fyud/Olk26U114qqjZtKufu7DCqn4tUZOe5BCYKTcnNqjRT18Ekv6dvLFjC7m84NH6/NHaO5eo3q9Ky5YeJKNDt3iHukc1K5bWa62W1N8trTDliAU8idPBzc+zn6sTmqHRVJTVXbysvKjUlpsO8nIeFSmf3wIwLBBUAAAAA3cug0ZyfMkadd5krdBxzMt9BKm2/zl1AoVCsWLEiJzvbWKca7Oq6ZfGjYkcBETlOi+YFNlnXqEqlPXersrxhNrZQwBvp59pLNozTGtSltTk3q86av8oTkRPPLZwfkzllhF6h5DhzxyTnCfr1s3CS5Fsf6wxqR55rVP9n2zpMozfkyVW5MqX51uBcDuPn6jTIw9m0Zq7WoL6t+CO3MpkNb3gcwUhpjJdzkLW3apHtgwpd1v7Pk03P1CVs7sLxTf7hXN4fn1JC7GYi/abHzQ6kjpcCAAAA9D03/rW+6KM1RDTogy98H3/S3s3pFXJyclYsX14tkxlUqgdHhjx/f7TYUcAIxQ7h4/hmm9xp9IbMwurrZfW/0HM5TKiv61CvXrFqVrOuCZaXc5ATzzWv+iwRuTsGBP7hmf3kTCJyHhUZ9uMx4rQ+Riur/AhbZYz3vI68/Zcp1bcq6vLkNTpD45u8gM8NcHfq5+roKRIQkUpXlVV+pLQ2p6Fhw4I9pnffaCgbBxXFqbsOZNYETImbGdxK0iBL+XrfFX7o/AWTJER06dDWEwVupqjDcikAAABAX1T5e8rlBVOISHzPzJH7frR3c3qF29evL16woFouJ6PxwVGh/5wzm+MucRg1hjdkKDH1+1gbiXJlNZeKqk0zlQdInEf2c7X7ntOtdk048lz7iUb2E49i39r/kB0sVF4kogEu47jvJsl2fUhEvqvWD/rHGy1PqNJVJd/6LxG5OwaM932sU40pqKy7WV5jvjk3ETk5cPu7OQ/1EjrxuaW1OX/IkkxdFoFuUwa6dssO7rYMKppFBUREJEtO+D6bRs5bNEnS5HNnSwEAAAD6ntrruZkx9+gqyjh8Cku56ThggL1bZGcGhUJ74fzZw788/+VXRPSPObNipk93GD2WN3CQ+WFlNZpz+ZXVdVo26erEjwhwc3d26O7m3Vb+kVvxu0pX2fEqXs5B/cQjm3UvaA3qM0W72KhjhMv9pbOXqPMuE0Mj9vzuNmlyszOcKdpVUZdPRFP6P9u1noQ6neF2pSpfXmsaIcbiMIybM9/ViavUZ9boLnG5VURGd8eAkdIYm3dZ2HLzu7JyaX//W4Kh5iGA1MPVaCyoqCByk2UVKBi30GCzYmnQAFFWZv412aQIslTasrPCfJAVIw41BR4WBl+x3SbsZ/cWkU+WgulMFQAAAIB2qIuLLy94QFdRRkT+b3xyl0cUxrq6ulOp+mvZRBTWv9+/FvzJxb//hD/N4/UPMD9MpdVfvN24+4QDjxPq6zrYs0NLslpDoSm5Un6EfbnviGZdE83wOYKx3vNSb2/XGdQ5yuOh8Z9ej4k0aulq3Lwxv2fx3NxMR5bW5rAXHeAyrssv+o48TqCnMNBTWKPR3SxX5VfU1qh1RGQwGuU15bwvjAAAIABJREFUGnkNEQUSBTKMnsOV16nKS6t/HCQZFCJtZ5/yTrFlUOEVHOwVHNw0T1ZexTBid3ciQ7lcSS5hTUIO8pSIiErk5WQgS6XUNKhgI4eGUVWy5ITv9ybQvEWTXJrkU3HqrgP7vqa5C8dL2cFU/afHPRZY36OyZ3vl1GdmBbGzOMr8p8UtCqo/c/tVAAAAANqhlcszH56uKcwlIu+lr/ovX2nvFtmZ9mo2G1EQEbf/wPvnzON6eTU7JqdUeamoWt8wT2CwpzDUz6W7N7PTGtS5FSnslAZqiBYsV3EReLc788GJ5zrGe/7Zol06gzrb/Y+Bqz8qWPO8trwk+9lnQnbvNR2WVX6U2IFJLVZ86gKhAy/EVxziK5bXaipqtBUqbWWtplJV3+FjNHL1OqleJ9VS8JVaulqc80BwgBPfNvPdbRlUtKTLSs1SMO6hwRKiYrmSqM0dScosljZ1+WByMdd/2sz6+EU6aYxPVnL+Ndkk4fUiRjxqfENc4zU0QJR5RV5OJJWVVzFc/6HstG+ONHJxXKRBViYnKjqUXEL9pzeECrzghyZf33rqbM74We5tVUFvBQAAAFimVyoz/zRTnXeZiDzmrwh8Z4O9W2R/WadPf7hz9yPTp0f/ZRVX0vx9SmcwZhRUmZZM9RQ6jA1wc3Hkd3erzKccENEAl3GB7pG2WoBV4hgwTDIjW35Ypassjg1w+XlWdeqhqt++L9zxP78lTxPRtYoUdqhVoNsU2676KnF2kDSMFtMbjJUqbUWtVl6rqajRKBpWjtLphApNVZ8IKi4fTC4m8h/b3qghhVzeXmnjGXRZV0uICRjcGB3ygh+KC5aVychz1soQoqZLSLE1pR6uxqyCY9sPkamrgSP1lNDllILGyIHl4SbWZ17NpdltVAEAAACwRK9WZz4Sq8o+S0TuDywa9uF/7d2iXkCv3/bdnnN5t5jzGQ+0iCjktZpTN+S1Gj0ROfK5o/u59nd36u4WqXRVf8iSTOOd3B0DhntE23wX6oGuEQpNSaHyYkVdvuvGZ/mx6drykvzVS90mR/IH9We7Rxx5rt00eZrF5TAeQgcPoQORkIh0BmNFrSZXfp3LVXsJw2x1le4LKtg3e9ewufXLwnpKRESqVg8VSySeZKm0A5eTekpNF2W4/tPiZgWxE8fZ4pBZK+nQ1hMFx+LjjxGR+Xq1+sZMEzfLVQAAAABaZdTprjw+v+b8MSJyuTcm+LMdba0ielfRFRX+np1DRMNCm4wsMhrpcrEiq6SaXTnI28Vx/AB3QTev76Q1qPOrzl6rTGGTPI4g2GNGu0OeuizYI7paU6LUlNyki0M3byx+8gmDlq699Bdm2/NsD8lIaUw3XbpVPA4jFQmkouE2Pq1tT8di3+aV5Go+6ZnjIRHRlbyr8vHSxiChTK4kcpF4EIcslXaMLDkhpVQ8cn4bq0WF1PdjsFMybh+NP0jLB5PFGditVEFcAQAAAG0xGLKWLan+PYmIhGOmDv9qD8Pr3qHmfcVviQfYD+OmTTNl1qh1J2/I2RH/PA4zqp9bD0zILq3NySo/alrfybbjnVrF5whGSWefKdqlM6hvDCvtt/xv8m2blad/rvungv4Z4+UcJHEMaP8svZ7tA0Fd1v7P9l1RUr/pcU22mOBIg/3Fxsq8LLOhTrKcPCUjDhgibafUHC94qDcZm46Yurw//tOfsuQVCnIbENxqfKDLupLbeIaH5oQKiRTlZcOH+De7KBlkKSez2q4i69TDAAAAgLtI9gvPVvycQEROw8aFfpvIFXTjq2rfknbqNBExfH74uHFszvWy2sNZpWxE4ebEnxHs1a0RhdagvlmVlnzrk/Mle9mIwt0xYFK/pcEe0d0aUbDEDt5sd4TOoC5fPkI0fgYROSal0o9/BHtEd/fVe4aNg4ri1F2fJ5dw/afGtfKLvvTe6BFCRebR1Pp3+EuH9mYpXEdHs90ElkvNjZgQKqrI3PNTFjUcmVzmP21msMRdTGYRguz3I5eVRFTfz5Dya0KqqcgUrgwfFyJUZO47lNNw8suJ+y7/kfx16tnEtqoAAAAAtJT7+mvlez4lIsGAEaHf/8QViezdot7CqFKdz8oiorGjR4vFYrXOkJJbfu5Whc5gZBgK9nGZPsxLKOiuLh2FpuQP2cHkWx+zE6aJiMcRhEpjxvs+ZvMZFBZ4OQcFut1LREpdacXbD+glAiISv5ugz77VY23oVrbc/K7pHhHmzDd/MJ9F3XKWguXS1q/F9Tet9Gpe3SVsbljRvuNl/tOemRXUdDMK8yrNapFpRVqLVQAAAAAa5W16//bG14iI7z1g1E+nBD4+9m5RTyuv0ZQo1Hwux5HHEfA5jjyOgMdlZ0dUnD9339w/EdGKVaseWvbs2bwKjc5ARE4O3HsGSjyE3bWl3W3lHy33vR7gMr6feGQP9E606nzJ3tLaHCKiK8XiuC2kJwffwWG/neO52ngrup5ny6ACAAAA4C50+4vtef9YRkQ8idfIg6edBg60d4t6js5gzJerrsmUpt2vm3Hkcw2/fLNmy7+JYf65/VvG3Y/NHygRhvV35XEYmzdJpavKqzp7W3nRtFAsEXk5Bw1wHWf32QvmO237/6SpWr+WiFzujQn99gD9P3t3Ht/EdS0O/MxotGu02ZZseQXbMhAwi4EkbDbZmwBpkjYvkLy2eW1M0ubXvgJt0qSvQJekrwW6vL6+xOmaNCZtkzbBQJImBJstBDCbwdiyDV6wbEm2bGm0azT6/TFGCC/Ci2xhc75/5KM765Xt8Jkz995ziPj/KCYSLh5CCCGEEBq9pu8/b/n9fwOAQC655e39N09E4fAGm7rcrXYPy8V6Q+0Lhm6TUsuM+XJDNh9RCAXkwixNujo+5RGYgCXI+b2swxd0BDmfM2CJropNkeJs5aKh6l5PvMiibVqkv+WZx+vOttnf/b3z0O7mn72U89yLie7dmOBIBUIIIYTQaATt9vqn/t15ZC8AkCKY9c8TygVFie7URGjt8TbZXN3uQGSLgCTS1dIsjVRAEgAQDoMvyPnYkC/Igb1r+mcfAsCF3IUtakMKLV6crZEKh1ny+BpBzt/O1PjYXmfAAgDRwcNAGklWOl04foli44Lz+8/cs9RrqgYCZry5T7vyjkT3aPQwqEAIIYQQGjHmzOm6L38+aGkBAHH2rILfvaWYfUM/v46d2882dXma7W5+RQQvSS7KSZJnaqRDTWQKnD4VOHEUAOSPf4WQjrKknZd1NPYcsnrqo2c0DUojyZJS6jzNshtkaOK6fG1tZ+6YEXJ5BArZ3E/qJJmZie7RKOH0J4QQQgihken8685L31kXDgIAqFY+MuN3rwtk415gIYEsjN9kdVmcvsgWMUVma+XTkmS05DoPk6H2NgAgtUmjiyiYgKXZccLsOhu9USPJAgClSE+RYqlQLaVUQlI8kamc4kiSmZlfVlH3+J0hl+fCEw/N/ddhcnJmIo7nSIXZbOY/GAyGSJMgiLS0tNjNcDjc0dER3yZJkqmpqbGbHMd1dnbGtykQCPR6fexmKBSyWCzxbVIUpdPpYjdZlrVarfFtCoXClJSU2M1gMGiz2eLbFIlEycnJsZuBQKCrqyu+TbFYnJSUFLvp9/u7u7vj25RIJFqtNnbT5/PZ7fYxNqVSqUajGVHT6/X29PSMsSmTydRqdeymx+Pp7e2Nb1Mul6tUqthNt9vtcDji21QoFEqlMnbT5XI5nc74Nmmapmk6dpNhGIZh4ttUKpUKhSJ20+l0ulyu+DZVKpVcLo/ddDgcbrc7vk21Wi2TyWI3e3t7PR5PfJsajUYqlcZu9vT0eL3eMTa1Wq1EIhlR0263+3y+MTaTkpLEYnHsZnd3t9/vj28zOTlZJBJd3atWt/7gBdubvwIAICBt439P2/jdrq6uQCAQOTguzZSUFKFQGLtps9mCwWB8mzqdjqKoSNMnVDZ1exzeq4uwU5USDcVqhKHIwVarlWXZwZsE8f73vnOorqG09Kn0++6PcbBerxcIBNFNR7D9gvUTV6gjcms5maoTzFWQqfzBFoslFApFzh1dMzU1lSTJ2M3Ozk6O4+LbTEtLIwgi0vT/6Xcd/7MZAKT3rNW/vI3f29HRwT+ox7HJP6XHbvZ7pB8mnP6EEEIIITQsfoul7suPus8cAACBQmb8/V7NiuJEd2pcsFz4Ype7weryBkP8Fookpicr8nXyES2HYFtb7rjnbsYXKH3mmWeee26YZw1MBWtQFE6iGU2jEQ6ffehzrs8+BICcn/7e8OX/SHSHRgynPyGEEEIIXZ/j0yP1X32Q7ekCAGnBoplvvDN5p7/H4A2GGqzuS93uYKhv4YRUKMjXKaYlyYSCERdNPr7vY8YXAICZVwppx9bsONHqPM6XqAMAihSnKwqzVYumcjjBI4iZfyg/XTI3aLvc8v2vKubOV86bn+g+jQwGFQghhBBC13H5lf9t+/Gz4RAAQNJDX8v/xW8m6cT3GJy+YJ3FdbnHy12ZxqKSCgt0dKZGOuoKCpX7KwGAEAqLrhdU2H2t52x7osOJbOWibNXiRFWpm3hCrbbg9V3nVy8Is1D/pTXzKs8ItdpEd2oEMKhACCGEEBpSyOczPbu+Z8/rAEAIIWvra+lPfi3RnYq/2k6mtsMZaaYqJfk6uZ4eUymJsNt1qr4eAIz5eTRND3VYkPOfs+2xeur5poRS5WmW3+CpYMeJct78nJf/dOk7XwnaLjdueHbmn8oT3aMRwKACIYQQQmhIF574gvPwHgAQJukL/rxbWbQw0T2Kv1NtjqYuF/85Wysr0CuUEuHYL9tTX9dgsQJAyV13D3VMs+NEU+8BPlEsRYpz1StyVFPwJzx8aU982XnsaPffX+n5cGfjc2l5/7090T0aLgwqEEIIIYQGd+lHm/mIQr5g5cw//VWUkpLoHsXfp5fs7b1eABAKyOL8ZLU0DuEEb//e9wGAIIlFxYMsZ7f7Wi90fxxZja2TGWck3TX1104MQ97Pf8laOx1V71pf30ElJeV894VE92hYMKhACCGEEBpE157dHb/9IQBIZ9025x/vT71FFMFQ+MjFbpvLDwBSoWBFXvJ1i06MSHtTIwAolKqiomsKjQc5f1PPwRbncb4poVRzUlZpJVlxvPWkJhCLZ/z5rZpVd3jOHTH/4kVKpc5Y//VEd+r6MKhACCGEEOrPfaG28eurAUCgVM16859TL6LwsdzBxi6+AIVSIlyRlyQZSa7Y6wrZrAsyDAf1KQ/+29ro7e2umrrujyKFsXPVy/I0y+N436lBIBbPfnvvmXuX+FtqW7d8g1IqU9c+kehOXQfWqUAIIYQQugbrdJ5aURi0tAABs94+pF6yNNE9ijNPgK1s6PYEWABIVoiXTk8SCkab4GkIgVPVgepjACB//Ct8Le0g5z9lebvH18ofoJFkzUlZhfOdYvB3dp697zb+7zD/1fdSVq9JdI9iGXG+YYQQQgihqYzj6p5cG7S0AEDm8zumXkTh8Ab31dv4iCJLI1uRF/+IAgBC7W0AQGqTIhHFsY43+YhCQqnm6x9ZnPY4RhSxiVNTZ79bSanUEIbGZx7sPXQw0T2KBYMKhBBCCKGrLv5os/PIXgBQrXw485vfTnR34szq8u832fwsBwAz9PTiHA056iIUMbBsyGoBAEF6X33Ac7Y9/Jpsg6JwSfpXdTJj/G86FUlzcma9fYAUQzgEdU+scJ46megeDQmDCoQQQgihPl17dne+8mMAEGcWzHjt9UR3J84u9/oONXazXJgAWJilmW1QjtON2Pa2t48ee/vYCUF6BgDUdX/Ml6HQSLLmpDxw89SziwvF7Dkz/nKAEALnhwuPLnfXXUh0jwaHQQVCCCGEEEDU4mxSSsx4c5dALk90j+LpUrfn6KVuLhwWkMTS3KScJNn43avuyKe/+te+X/1r3xlzR7urhk/0pBDp5+u/MH43ncLUy5YbX3sfCAi5POcfKfE2Nye6R4PAoAIhhBBCCFins3bt/VwAACD/lffl+VNkfo6f5Rpt7qrGrurWHrhSjCJVOaZS2ddVuf8TACCEwpTp9DnbbgCgSPEC/SM4RjFqSffeN337GwDA2q3nPl/i7+xMdI/6w5SyCCGEELrpRS3OTn3mv5LuuTfRHRqrQIhr6/Fe7vXaGH9ko1QkKM5LVojH9/Ev7HYdOHMWAHJzcxrc7wMARYpxWfbYpa59wt/R0f7z7wYtLec+f2fh3oNCrTbRnboKgwqEEEII3ewubv0vfnG2YtHd07+/JdHdGb1giLvc67vc67U4fdHbk+WiDI0sUyMVU+M+S8V+obbBYgWA9EIlX49iRtLdtEg/3ve9GWRv+A7bZbP88ef+ltozdxbNfOt9ecGMRHeqDwYVCCGEELqpde3Z3Vn2EgAIk9Nn/fmvQE7KyeFtPb7WHneH45pYQisXZailmRqpNK6F7WKrrqwCACAgd1kaABRo705XzJmwu095uS/9LGjvtr/3h0Bnc839s/Nf2Zt09z2J7hQABhUIIYQQuplFFmcTQpjxl92URpPoHo1YMMQdvmjvcl2d5qSWCTPUsiyNRCZKwJNe5aGDAMBSYeOcdIOiMEe1cOL7MLXNeOX3zXlG8/bnOU+o/kv3Zjy3Les/Nya6UxhUIIQQQugm1rjh6/zi7Gk/fZ2eOy/R3Rkxl5891NTt8rMAoJIK09XSLI10vFdNxBCydJ6+1AwA8xZNU4j0c1IeSFRPpracTc+pFt/WUPp51tF7+b83eS6cz3lxiyQrK4FdmpQDfAghhBBCY9f6y+3u01UAoH/yO6nr/j3R3RmxbndgX72NjyhmpirvnqGblUonMKIAgJ6L1fcVFirEoiWrFy9OezyBPZnyNCuK5+6vlc68FQDsu/545s6Z3R9/lMD+EOFwOIG3RwghhBBKCObM6XMPzA+HQLH4nsJ3P4DxKCw9npq7PdVtPeEwkASxOEeboR7fLLHDwQQszLt/VzoFHnlY+sgjuDh7AnCBwMXvP2994xcAAASkPbt12vPfT8i6IAwqEEIIIXTTCXk8p4rnBy6bBHLhvMNtYv1kevwNA5xtdzRYXQAgoshluUlamSiB/QlyfiZg6fG2ttmPLdwvBYDwTCO99M4EdulmY/nbW5eeW8v5AAAUi++Z9ae3Jn51EK6pQAghhNBNp+l7mwKXTQCQ89O/TK6IguXCRy52Wxk/ACglwmW5STLRxGV24nlZBxOwMn5Lt6/Fxzq9bC+/XdvVF9tIs6dI6cDJQv/oY3TRwguPP+hvqXUd+9ep4sKCP72rXFA0kX3AkQqEEEII3Vy6P3i//sn7AUDzucdn/uEvie7OCHgCoYNNXYyPBYBUpeT2aVoBGbdZW0zAEuSuppDysg5f0BF9QJDzOQOWHl/rUFdIOyWsevfU/GnTlv1gKwgmOtRBIbe7/htP9X64k2/m/mqn/tHHJuzuGFQghBBC6Cbit1hOL88LMS5hctr8I/UUTSe6R8Nl9wQONXUHWA4AjDq6MF0Zx4vX2PaYXWdHepZCpJdRKlqk00qzaZH+L9/c8Ovde2il8sDZEV8Kxcvlsv9r3fx1/rNy+erpL++Q5eZNwH1x+hNCCCGEbhrhcH3pl0KMCwDyX/nbJIooLvf6jjXbuXCYIGBRtjZLI43jxYcfUWgkWUqRXkKplWKdVnJNAlPO4Th49iwApKWnx7FvaKQySp+hFyys//f72F6782DF6WUVqU9/P/s73xPIZON6XwwqEEIIIXSzaPu/37iO/QsA9F99Tr10WaK7M1xNXe5Tbb0AIBSQS6drkxXiOF48ElEoRPqZSXdFtgtJ8YjSN/WY6k61tgHAwiVL4tg9NAqqhYvm7qu5tPkF++4/A0DnKz/ufvuPWT/Yrv/iv43fTXH6E0IIIYRuCu7a82fvmx0OgjjnlvlVJ0lRIjMmDZ+f5d4/38lyYbmIWpGfJI9fkewg5z/W8aYrYAEAMqhq3A8CgioqKioqurrAt6qqqr6+vt+JNE2vWrWKvjLOYzabd73xetvhw3vP1hAk8drf346+AkogZ/WJi89903P+U74pK1ye/6tX5TNmjse9MKhACCGE0NQX8vtPlxT5m88TQij84Jx81i2J7tFwnTU7TRYGAG6frk1XxW3WU3REoRDpf/OdfzWYGgGApun9+/fzx5jN5jVr1gx6+oYNG9atW8d//sG3vlXx3nv8Z4FKeeIMLqi4kYTDlr//tfUnzwWtrQBAkJDy+Ldzvr+FUsZzTQ5gRW2EEEII3Qwu/df3/M3nASDju7+cRBGFn+WabC4AUEmF4xRR6GTGD1+7xEcUALBgwYLIYQaDwWgcJDmsQqEoKCjgP7ONDbPDQYVYBAQIVKrSbzwbr06i+CAI/aOPLThmMmx4iRRDmAPrG784uTjb/Oc/QFyHFnCkAiGEEEJTnH3/J3Xr7gQA+dwVc9+vnETFs8+2O01WBgCW5ialKeNTMzs6ojAoCpuPclu3bgUAo9FYXl4+oksFas8HjhwAAKCE0s+tFkyqih83IX9nZ/OPftD9j9f4pnzBytnv7BVI4vN3hUEFQgghhKayQFfX6RUz2Z4ugVw4r+qSePLkJvKz3N7znSEurJWJ7ihIics1+0UU4p58fhaTQqGoqKigR5IOK3DubODoYQAAoUh63yqMKCYL5+lTl17c6D65HwDmHqiX58enUiFmf0IIIYTQVGZ65km2pwsAcl56fRJFFABQZ2FCXBgAZhviM/2dCViOdbzJcn4AMCgK56Q8UPpiKb+rrKxsRBGF/8Sx4OlqAACRWPq51YKU+MQ8aAIo582fu+cT63vvBrus8YooAIMKhBBCCE1hl176ofPQbgDQ3Lt2IqsLj50vGGqyuQEgWSHW0XHIITswogCAkpKSkydPbt68OcfNuP70O1KuEOTmC/PyyZireP1HDgVrawAARGLpA2sESclj7x6aYLoHPx/fC+L0J4QQQghNTZ07/3Jxw78DgDAlY/7h2klU6g4ATl92NNpcAHBHQYpWNtbst+2umrruj/iIIle9LE+zPHpv4HxN4NND0VvIpGRqWq4wz0goFP0u5T9QFTTVAgAhkUruXyPQasfYNzQ1YFCBEEIIoSnIXlVZ99hKACDFcMt7p+i58xLdoxHwBUN7z1u4cFhHS1bkJY3lUu2umqaeQ162l2/OTlmVrpgTfQDbaPJV7hvqdDJF1xddyGQA4N33cehSA/ARxQMPCjSasfQNTSUYVCCEEEJoqnGb6s/dPyfkDhIkFJTv1xaXJLpHI3OqzdHUNaZhiiDnb3Ucb3Ye40cnAIAixTOS7k5XzDGZTGVlZaWlpUajkW1r8X24FwBAJJau+rxAqw17PMGmBrapgeuyRV9QkJoWdrs5xgkAhEwufeBBUqUa49dEUwmuqUAIIYTQlOK3Wmu/eE/IHQSAadv/MukiCk8gdKnbDQCpSskoIgov62hxHG93nY2EExJKladZrpMZhaSYYZhNmzaZzWan0/l/P/mJ7+N/AQAIKOm9D/ATmQiZTDRnrmjO3FBvL2uqY5sawm4XAIQ6O/irEXKFdNWDJB3n0mlossOgAiGEEEJTR8jnq/3i/XzxYMOGl1IfezzRPRqxOgvDhcMAMGeESZ+8rKOx55DZdbWgtUKkz1Etip7vxEcUAFC8aJH3/QoIsUCSkrvvG5gQVqBWCxbfJl50K2s2s40mtvkiBAOEgpat+vzAhRYIYVCBEEIIoamC4+r+4wmvqRoAkh4pzfnO9xLdoRGLDFOkq6UqqXCYZzEBy4Xuj3t8rZEtGklWnma5VpIVfdj27durq6sB4IF77vm8XBz2uAFAcsfdVEbmkJcmCCo9nUpPh6XLQzYbqdUS4jikokJTDwYVCCGEEJoimr7/vGP/OwCgXHK/8Ve/TXR3RqO2k+GXu85OG+4whd3XesrydmSyk0FRmE7P6RdOAEBFRcXOnTsBwDh9+tcL8viIQrziDipn+rBuQ1GCtLRhdgndhDCoQAghhNBU0P7H31n++HMAkOQvmPnG3wmBINE9GjFPINRidwNApkZGS4b1kNbuqjln281/NigK8zTLpNQg66dNJtOOHTsAQCGX//ieO+RBPwCIblsmNBbErffo5oZBBUIIIYQmve4PP2h58SkAEKVNn/3OBwKZLNE9Go3zHc5wGAiA2WnDKqlRY9vDr6CIZHYa6sjy8nKGYQDgJw+u0UEYAITzF4pmD3k8QiOFQQVCCCGEJjfX+XMNpZ+DMAhoxay3/yVKSUl0j0bD7Wdb7R4AyNLK5OLrPKEFOf8py9v8IgqKFC9Oe5wW9V9pHW3dunVms/neDMNcpQwAhDNmi4sWxa/vCGGdCoQQQghNZn6z+ew9C4LdFkIIs/5xTLVwsj4rH2vpabV7CAI+NytVJoo1d4sJWM7a9rgCFgBQiPSL0x4XktdfPO39ZF/oogkABNON0jvujFe3EeLhSAVCCCGEJiuWYc5/8b5gtwUIyP/trskbUbj9bJvdAwA5WnnsiCJ6WbZBUTgj6a6hIgqGYRoaGhYsWBAOBLwf7uEsnQAgyM7FiAKNBwwqEEIIITQp9R46ePH5b/ou1gBA1g9+k7xqdaJ7NHrnOpgwAEkQM1NjraaIXpadq16Wp1k+1JEmk2n9+vUMw3x7/fqHlDK+ErYgLV169z3x7TlCPAwqEEIIITTJ+M3miy9s6vlwJ9/UfWlDxtPfSGyXxoLxBS/3eABgWnKsYYrhL8uORBRhv1/bVM/l5QKAID1TcidGFGi8YFCBEEIIoUmD8/vb/vfXHb/+Ll+VQZicnv3DX+seejjR/RqT8x0Mv8J1pn7wStXRte2uuyybYZi+iMLtvjc/d3leLgAI5y8SFy0cj84jxMOgAiGEEEKTQ/cH7zf/4Nv+tnoAIEjQfeW7OS/8QCCXJ7pfY+Lys5d7vQCQl6KQCK8ZpmAClnamxupp8LK9/JbrLsvuiygcDo5x3jdzxgtrHgCR+DqkaxmnAAAgAElEQVQ1sxGKBwwqEEIIIXSj8zY3N238hvPIXr4pnXVb/q9fU9wyO7G9ios6iwsACAKMur5hCruv1eo2RccSvNjLsuFKRGGqreWczvtmz3phzQOkSiO5935SOdzi3AiNGgYVCCGEELpxhbzelp+/bPndj8JBAAABrch88VeGLz0JBJHorsWBN3i1hLaLvdjkMFk99XxmpwiFSJ+uKNTL8wctlR2tsrKyvqaGY5jPFc5+Yc0Dgsxs6Z33AIUPe2gi4N8ZQgghhG5Q3fs+vvitJ4LdFr6Z9Ejp9B++LNRqE9urOKq3uPiCYY7Q26csXdG7dDKjRpI9nFiiD8fdLhHdOy07VaX+j5JloqJbRfPmj0OXERocFr9DCCGE0A2ns/yN9l//1N9SyzfF2bNyd5SplyxNbK/iKxDi9pzrDHFhocgskX3Mb9TJjDq5USczDqeeXTRPxT/5ShQgEkvuvJdKT497hxGKAUcqEEIIIXSjCPn9nW++3vGbnwY6LvJbKE1y6vrvZX1rQ2I7Nh5MVneICwOAUHwaAHLVy7JVi0caSzAM88nf/35rOKgIcwBAapKk93yOoGMVu0BoPGBQgRBCCKHEC7lc5j+81vHqT1m7ld8i1GenfeP5tC89KRCP7Dl7UmC5cKOVAQAB1SkQdGkkWTEq2Q3F2dX11L/9m6mpaZkx/+VHHxZk5UjvuBsXUaCEwD87hBBCCCUS29PT9spvrH/4acjl4beIs2elf/N5/aNrian7fNxkc7NcGABE4nMAMIqIoufsmfVPr28wdwKAQaMWLlgknrcASDLuXUVoOKbs/6sIIYQQusGxDkfzyz/qKt/OBfu2SAsWpW94UbfmwYT2a9xx4XC9lQEAUmCnhGaDolAryRr+6WGv1/avD5596eUGixUAHrj91udeKSNVw1vPjdD4wKACIYQQQgngt1jO3r0waLvMN+XzijM2vph0192J7dXEuNTtCbAcAIgkNQCQp1k2/HPZhvquyk+++bs/NlisQBCrHnjgR7/5zXh1FKFhG6+govPIm7vOCeY99NjilH57at8rO2gBPrd0+p2lD+SOYC9CCCGEpoKQ21376AN8RKEqeSj92Q3qpSN4sJ7UwgD1FgYASJIRClty1cuGmTQ27PH4DuwPXW7lIwpCJFrz6KNbfvzjce4vQsMS/6CCsx18658XXAAA/cs38ruEs79YukQLAOf3vrKv7K3uK4FH7L0IIYQQmhrCweD5tQ95TdUAkPnCLzL/338mukcTqq3H6wmEAEAkqaFIcbZq8XDOCtRdCBz7FAL+Bou10dZFKpWrH354y5Yt49xZhIYr7qt5aiv+eTnvodL/WKEfsMt26ONaNz37ziV9NWtuuf+RGbTjzMdH7NffixBCCKEpou6Zr7mOfwQASV98+maLKACgtsMJAATpFYqaZiTdPZwcsv5PDwcOVULADwCzSu7YUVa2/tlnMaJAN5S4j1TMerB0FgCw3f13cLa6ywyhnj0jqgxmijFbUXeutdG2ZCHE2jtwsIKte+8PB/rqaxL07EfW9sUi0dsBlNHzrzqPvLnrnJv/rJn9hS8uidzKdmDnP+oYYiSnIIQQQmg0Lr30w549rwOAcukDxl/cdIsBzA6vy88CgEh8TkIp0xVzrnsKe7kteP4sAJBqjXjFHQKdrgSgZNx7itDITNxCba7b7gLlvPxrnsuTtQoAi70bOIi1F64NKvjIIWtF6X0zgA8J3tkJj6xdorxmO3QeeXPXP9+Chx5bnALn975y+HLmnaWP5/bNs3r7970lX73fyK/i6MpYWbrW2Hfl65+CEEIIodHofOvNjv/ZDACS6XNmvv43QiBIdI8m2oVOFwAQhF8oMs1Jeey6x4cDAX/VJwDw9vGTt61/ZoZON+5dRGhUJi6o6LK7AIb8tyP23mvV7jnQKchYyUcOAClL5qfWHWhttC2RX+wg6MLFfdtBl5+lOHfB3g2QYut2EIKMfH7ZN5myfF3pcs7WZQfo2HvAApl3XgkVqBkPLr34ytHjpsX3a4Y6BUcrEEIIoVGwV1Ve2vgEAAiT9Lf8/QOBTJboHk00m8vf4wkAgFB8QStNH04aWd/hg3s/++yPVYcsIe7Ppsb9+/ePfzcRGo0bIqUsY4+1boKxX/Mkz9Y1WIDImn51xICa8WDpDFuXDZLvf/oWgGtTSPFnpiSpwnWXK3+/FyJDDWRKshZqD16+GjnwktR06FxDEzwwxCkIIYQQGjHmzGnTV1aGOSAlMPOvH4sNhkT3KAGuDFOwInHtnJTS6x7/3h/+8OqvftnhcBJCIalSpaWljX8fERqliQsqkrUKAO+gu2itNhli7R3G5VOSUyASTggyVpbeb+TTSfG7b7n/adj7yuHLlWVllQAQna82dHVjhDr2KQghhBAaPl9b24XH7uR8AAQY//iJ4pbZie5RAvR4AlbGBwBCUV2Oan7sNLKVlZU7tm1rO38ewmEgiPT8/NJnnlm9evVEdRahEZu4oIJM0irgQkuDfXHK1SChy+4CUGqTgIRYe4fHdmDnQSs95wtrB19PfUvfOAa/JKN9X9keeGo6xFyBPcgpGFcghBBCI8H29p7/wj1srx0Asn/ymrZkZaJ7lBh1FhcAAHAyaVOu5quxD960aRPndEI4nKai1z/7zQefemoCeojQWMQ9pezQd0qZkUGHe1vqoqY62UwtLoLOyku5zt5o1Ix8PYSvnTFV+17Zqx/U2XsYUGfPGDQ+YOsuNF29woNrZssBmO6umXkZ/W4KnO3gp3VDn2Ib1ZdHCCGEbkpcIHD+sQf9rXUAoP/a8+lPfi3RPUoMl59t7/UCgFDcmKOec900svcvWZIqFX9vzf3v/PznGFGgSWHiggqAlGV3zZIz5/ZdqTxxfu87dYxq7l38MEHsvdFm3Tpb0XPu7Q/q4MqRB7oyVt43Q6uhISpCsB36uJZ/J8DWvfeHAwc/2RkpeXE1XJm56BY5c+6fe01XLl5b8c/amgNvHTleMdQpCCGEEBqWcLjuqS+7zxwAANXKR3J/+FKiO5QwtZ0OAAAI07KWPM3yfntNJtP69esrKir4Ztjr/e7CeX/7f19/YPFiSfFNOrCDJh0iHA7H9YLXLJLmCTKiM7FGHzBwlULsvVdF16OIun706cp5D83r+GdVV8bKr95vvLYYRYwuQSQjbcxTEEIIITS4kNvtOnuGOXuaOXq454NyAJDNXjKnYp9AIkl010ap2x0IA5AECAiCJIAgCJIgSAJIkiAJgiKJ2Kd7g6E95zoBgBJdWjJNp5Nd8zjBRxQMw6SlpfFxhfeDPaHLrQAguW8VlZE5bl8LoXiKe1CBEEIIoZsOc+Y0c+aUp+as49BH/ubz0btEhtzCD4+KkpMT1bex8ARChy92O7zBuFwtLfn40szPR2+JRBQAsHnz5tWrVwfqLwQOVgKAcMZs8bL+YxoI3bBuiJSyCCGEEJp0LH97y3Wqmjn5mefswYF7STHICu9WLFic+qUnJ2lEYWF8Ry/1BENcXK5GiS7NTFkQvWVgRBF2uQKfHgYAUqkS374kLvdFaGJgUIEQQgihkel8q7x9xw/9bfXRG0kpIZt9l8iQSS+6lS5apJw3P1Hdi4vaTqa2w8l/zkmSS4XXKdHLz4YKh0Nu1upl7YzfHAr7geAAOAI4IIK52uzoandms7lfRAEA3v0fAxsEghDfcQ/cfOXG0aSGQQVCCCGEhsvyj3cu/3xzZIKTMEmvXf3v8rnz6Xnz5TNmJrZv8cJy4c+a7R0OHwBQJLE4R2NQSWOfwgQs3d42q6e+x9fKbyGFfclwFCK9XpavlxtpkT76lE2bNvERxYYNG/iIInD+HGfpAABh4XzB5BzbQTezeAYVZrOZ/2AwGCJNgiD4ApAxmuFwuKOjI75NkiRTU1NjNzmO6+zsjG9TIBDo9frYzVAoZLFY4tukKEqn08VusixrtVrj2xQKhSkpKbGbwWDQZrPFtykSiZKTk2M3A4FAV1dXfJtisTgpKSl20+/3d3d3x7cpkUi0Wm3sps/ns9vtY2xKpVKNRjOiptfr7enpGWNTJpOp1erYTY/H09vbG9+mXC5XqVSxm2632+FwxLepUCiUSmXspsvlcjqd8W3SNE3TdOwmwzD8g04cm0qlUqFQxG46nU6XyxXfpkqlksvlsZsOh8Ptdse3qVarZTJZ7GZvb6/H44lvU6PRSKXS2M2enh6v1zuaplh88S+vd/3vz9jWvqqyApVG8aVvq9Y+kZyWJpFIenp6HGYzAGi1Wr7Jn8s37Xa7z+cbYzMpKUksFsdudnd3+/3+sTR9Iahzhj0BAAChwKumq2t722vsfhgJkhApiFSlIGta8jyZhO7q6mK6AgyYk5OTRSJRV1dXIBDwer3hcHj16tX33HOP2WwmPR7FsU8BgFOpuw0ZYDanpKQIhUL+YAAYtGmz2YLBYHybOp2OoqjYTavVyrJsfJt6vV4gEMRuWiyWUCgU32ZqaipJkrGbnZ2dHMfFt5mWlkYQROxmR0cHv/g5jk3+KT12s98j/TDhQm2EEEIIDY3jLP/8R/v2rb5L5/gNlCY57evfN3xt/eTN5hQtyPnbmZoQ5+v2tTg8MiczPxwWAgAlapZIDxNEaPiXklAqvcyolWb1y+80KIZhOjo6jMa+Iz3vvcPZrECS0of/TaBWj+67IJRAOP0JIYQQQoPhOMs/3m7f8aOr4YRWl/bMC1MmnPCyjsaeQ2bXWb7p984L+AsBAIATS4+LxPUUKaZF6UqRnrperTqpUK2VZEkp1fDvzg/o8Z99VZWczQoAoqJbMaJAkxQGFQghhBDqz/HpkYZvPhm43FccltLq0r7+ouGrpVMjnGAClmbHiUg4EQ6LvJ4VoaABAEQUO8PgTpYv1Eoejvt9d+zYAQAbNmyI3hi4UMs2XAAAUqcXzZ0X95siNDEwqEAIIYTQVWxvb9N/Pd/99qt8c4qFE3Zfa2PPwchyagBQC+d1OQpDQQAAvVJya45GJCDH49Zbt27la9sVFxcXFRXxG0OdnYEjBwAAKKF4ecl43BehiYFBBUIIIYT6dJa/0frDb7KOXr6pe+I/p/3o5akRTlg9Dc2OY9HhhEFRmCS+/dglH8uFAWBWmnJWKj3ouSdPnszPz4/MVuKZTCZ+0X+EQqGIrJHgMQzT0NAAABUVFXxEYTQaIxFF2O3yfvQ+hMMAICm5U6DRjv1rIpQoGFQghBBCCLzNzY3ffoY5+gHflOQvyPvVa8r5C2KfNSm0u2qaeg552b5IiSLFOllBnmaZPyCrbLDxEcWS6dpB88bu3r27rKyMT4aza9euSDKcjRs3VlVVDTx+1apVW7Zs4T+bTKZ169ZF7zUaja++2jcEBCzr/WAv+H0AIJy/iMqZFo/vilDCYFCBEEII3dTCwWDb//7a/ItNXAAAgBSBYcP2zGe/RUz+4mt2X+s5257ocCJbuShbtVhIihlfsKqxi48oFmcPElFUV1eXlZVVV1dHtkSPS/AZ6gfiEyXzInk5eXxEERnu8O7/mOvpBgBB9jRx0cJRf0eEbhCYUhYhhBC6eTlPnWx89iu+izV8k77tvrxf/J80JyexvRo7L+uo695n9fTV/JZQqnTFHD6cAACXn/3EZAuwHAAUZWmmJcmizzWZTNu3b4+EEwqFYt26dSUlJdFTmxiGMZlMA+9rNBqjZ0mZTKZImBG9y3+yOnjyGACQmiTZgw8DhS950aSHQQVCCCF0Mwq5XJd++APrX34BYQAASpOcveV/9I8+luh+jVWQ87c6jjf2HuSbFCnOVa/IUV0dCvAEQp+YbL5gCADmZajzUuT9rlBaWnry5Em4Ek6sXbu232qKMWLbWn0f7gEAEEvkD3+RkCvieHGEEgUjY4QQQuim01n+hvm3231NZ/hm0hfW5/7wZUqjSWyvxq7f8gmDonBG0l3CqCoT3mCoqqGLjyhmpSkHRhQAsHr1an45RNzDCQDgHL2+Tz4CACAI6d2fw4gCTRk4UoEQQgjdRLo//KDlJy/6Gk7yTXHWjNwdr6mXLktsr8aOCVgudH8cSe6kkWTlaZZrJVnRx/hZbr/J5vKzAGDU0YXpSgAwm82vvfYaTdP9ykeMh7Df7333HY5xAIDo9uWiW2aP9x0RmjA4UoEQQgjdFFw1Zy9tfp759H2+KVCq0je+nFH6TGJ7NXZBzl/X/XGkkh1Fimck3Z2umNP/sBBX2dAXUUxLkhemK6urq3fv3s1negWAVatW9UsIG2fhsO+j9/mIgsqfgREFmmIwqEAIIYSmOF9ra/OPfmDf82d++QQpAn3plqxvbRQoJvfcGy/rMDM1zc5jLOfnt+Sql0VWY0djufCBxm7GxwJAllbWcfbQ+p9URGd2Ki4uHt+IAsB/9EioswMAyBS9ZHnxuN4LoYmHQQVCCCE0ZQXt9tZtP7W+8fMwCwBAkJD0hWdyvr9VlJKS6K6NXpDzWz2mduZsdCU7jSRrTsoqKaUaeHyICx9s7OrxBAAgQyP75fe+EZ24adWqVaWlpZECFOPVZ1N98PxZACBkcum99wM5LkW7EUogDCoQQgih6+NLIy9YMGmKwYV8vsu//Z/O374Ycgf5LaqVD+dseUluLEhsx8bC7mttZ2qsnvrI0AQAKET6fM0ynWzwcQYuHD7Y1N3tDgBAqlJya7aGLx8xTpmdBsVeuug/8AkAACWU3Hs/MSUqlCPUDy7URgghNBFMJlN07TBefn5+9COd2Wzu7Ozsd4xCoeg3L4VhmPF+EDSZTGVlZXyFgehJMkVFRZGKyAzD7Nixg6bpcZ+LP1LhcOff3mp76fmgte9Fvmz2kmk/3q669bbE9mvUvKzD4m5odR6PpHUCAIoUpysK0+k5tEg/1IlOX/CPuyov1JyeVXTb/DmzlucmAYDZbK6url69evX4dppl2daWYEtzqLUZggF+m2Tl3VRu3vjeF6EEwZEKhBBC11dVVVVZWVldXW02m1evXr158+bIrq1bt/arHAwAfC6dyJSS7du379y5c+BljUZjeXk5/9lsNq9Zs2bQu2/btq2kpCT6UqtXr96wYcP4hRbbtm3jKxX0E10vuaKigl/jW15ebjAYSkpKboTowll94uJz3/Sc/5RvijMLMr/3ku6hhxPbq9HhpzlZ3Q2RGnY8ncyYTs8ZamiCZzKZ/vTXdz76+BO7zQIAl00131jzZ36XwWAYv8lOYb+fbWlmLzWF2lqitxMyufCWQowo0BSGQQVCCKEh8elxKisrox+m9+/fHwkqqqurI8lz+ikuLo48ukWfPkaVlZUAUFFRUVlZuXHjxlWrVsXlsiaTSaFQRDq8bt26yK6FC/vqphUVFUXHDCUlJZWVlXzsYTaby8vLy8vLaZpeu3ZtaWlpXHo1Ir62tkubv9fz/pt8U6BUpW94KWP91ye+J2PEBCxWd0Onx+QKWKK3K0R6fmhi4Drsq+cyzO7du9/bVXH6XG0wxPEbJULBugfvE5DE+PU57PEEL10MNV8MdZohag4IIZVROdOp6XmCtLTxuztCNwKc/oQQQmhw5eXlO3bsiN5SXFxcUFBQVFRUVFQU2bh9+/b6+vp+5xoMhujFrwzDRC+NjYi+zlCHpaWlRb9XNpvNW7ZsiQwjFBUVbd68eSwvnqurq8vKyqqrqw0Gw65du0Z6utlsrqysjEQXvBMnTkT2Pv300wNHcgwGwyuvvBLpdllZWVlZWb9jaJp+6qmnosOboYQ8ntYdP7OUbeWCAHxyp6c2Z/3npsmV3MnqabC4TT2+1ug5TgBAkWKdrCBHtTDGNKeIsrKy377yqtvPcuEwAGROy336yS89cO9dYxzUCnV2hr1uzu8Hvz/sD4Dfx/kDEPBBuC9uCXV1ReY4AQCIxFT2NCrPSKWnj+W+CE0iOFKBEEI3hOrq6ui5+zyDwVBcXBx5HuLfwg58619SUhL9Br2qqmrQp/zoS5nN5qqqqtiXijzyFhcXl5SUlJSUDPpktnHjxut+O5qm+8UPoz7MYDCUlZVVVlZu2bLF5XJVV1c//vjjGzZsGMUU+Ug4wTdH95bNYDCsW7du3bp1DMPw0UX0V+jo6BgYUQCA2Wzu6OiI/IQjQUg0/td9naCC4zrK32h7+bus3QoABAnaR57OeWGzODV1FN+FF+T8TMDi9FtDnO/wsSqKFMlFSZG9OdMzVEqtUqzjm40NF48c/KzfFVJT9UtW3KpQ9BWrdrncRw581tlp6XfYkuW35uVP97KOgROcAEAh0idJsmIvmWAY5uTJk/n5+fxPMhwGUVIm4wuqk3XLP/fQXXeuvG/hTGpsAxShzk7/pwe57q5hHU0JqWl51PTpVGbW9Q9GaGrBkQqEEEo8hmFWrlw56K4NGzZEHiuHWplA0/T+/fv5zzFWJkRfasuWLbt37459Kf5q451qc9QYhikrK+N/IAqFgp8WNUwmk2n79u2RcGJcEwFVVFR0dHT025iWlhYdBZnN5oG/Dn75SiREqays3LlzZ1FR0apVq/hfSu/hQxef+3++ptP8ASNK7hRZE//Z8aP+EONlnSznLSrJFWv8fGKlo/suvPnrTwaeKJWLf1b+Nf6zx+1/bt3vBr3+w19dvnJNIf/5ndcOVe4+E/tSAHD5YhcALJi9RCPJ1svzB00Oy2MYhl/kw//S+T9aTyD0WbOdz/JEkcS8DHVOkmw4P4qhhD0e37HPQo11A3cREikhFoNYTIjEIJKQYjGIxYKUZCp72ljuiNCkhiMVCCGUeDRNp6WlDXz0jJ7lDwBDLQLOz8+PfDYYDENdqqDg6uNmUVHRoEFFv5SpN2xEAQA0TW/cuLGkpKS8vDyyjHtQ/FhB9HdZv349P0ozAXlFhzOEws8Wi31MeXn58cOHj3700f9s2ZxFkYud5vkiSBYAAMhuuT1n68/US5fBYBm0Tpw4ET1jjQlYvvTEfzSYGkPhABdmw3D13WKbde4jTy3jP7df7AYAASEMQ5jji1wAwLXjOTK5WJtC222DLJjR6q7OvEqfnjTwAADQpvT9zClS7O1Q/Po75RQpSjc0lJSUFBW5i4sHKQ/Hr/DpF0CmpaW1O7wnWnr5RRQqqXDJNK1cPKYnnMDpk4HTJ4Hty8YrnDGbmjmLDx4IoXAsV0ZoqsKRCoQQQlNKZWUln762vr4+evbRq6++GnmwXrdunclkmpiqZ6Pmt1jc58+5a8+5z5311p+rrzn2exfZxl5zDCkSffXLX974f6/wzYqKiq1bt/a7TijMhsLB/33jB2Ity1eL+87a13yeQL/DtCn0M999fMGC+VJKqZVmh7yCo4dODvzhGI1GSsp5WUe8vqaQFNMivclk6jfRi6ZpPqFW5LfW79spFIqSkpIVxcXqvPkXu9z8xnydYo5BSRKjn/LEtrcHDlVxTN8XJHV68dJiQdLgcRFCKAJHKhBCKAEYhtm6dSvDMNu2bZuA2lsD8fPm+c9e1uELDvcZUUBKIvPptZIbbuK4yWTatGnToLuiF5CUl5dPQLGLUet6f6/t7+WRPE68LCFs1XBdITjsg08FdLdASCloYXIyo0+N/DZdQZs/1FcMhJ/FFAoHAUCbQntkFz2+vkt96ycPXb5ky0yfJhdqZhTMzEyZKaVU/acbSWD16iEXGceYmzQ6RqNx165d0UveGYbhk/bu2rWLj2346IKPJXgtdm+TzcVHFCKKXJytSVWOvq5c2OXyHT4Yamvmm4RMLl50K5U/iWsFIjSRcKQCIYQmWnV19aZNm/hn3OjX5+Oq3VUz6HLYuFCI9HyWTymlllEqiVAlpVQJCTnMZvO6detcLhdfMo+m6YKCAoVCsXDhwoRXkLiu7o8/6nrvnZ69v+M8oejtArlQUrBMNHsmlTedMBrAqPdS3rr6utozJhB7C2+bLpNfTbFqqmnvd9mM6cn8ARpJllKkp8V6pUg3nExKiRJZ8s6n5IpUMongwuHmbm+9lXH7+wZuUhTiW6dpJJRglLdkWf+pk8FzZyDEAgCQpHBWoXjhIqDw3StCw4VBBUJovPC5Qfs9yUWXVU5NTb1hZ56Mn507d27fvp3/vGrVqi1btozr7byso8VxvN11ln9vPcEoUkyL9EqRXkKplWLdIK/Db3pBzm/e9y6z91+uir9w7qtTkojkJPGaB4Oz08J5OreOHcWvj4/0JkUUMXwhLtzU5TZZXb5gX9wlosgCPV2gG33yXPZik//o4bCnbwKVwJAhXrqCVOEfKkIjg0EFQig++IyoJpOJLzUQmWoSnXFo0GT8AFBSUrJt27aJ62uC8FOeImtMN2/ePIosqMNn9TS0MzXRQxMUKU5XFFKkGAAoUhqZxTSiZ30v6+Dn07Ocz+nvm0Bl97XyH5iA5bqPv5oBIxj8+EakmcCBjpGK/DR4Pd7WyGcP6+hXbwEAmIAl2NMLTh/h8ILLK9xXK6r6lHRdPSCkFftX3ha66xaYmxHjvhJKJaVUQlJCi65MRZNm932YDD+3UQiEuEabu9HmCrB9pSGkIkGBjp6WJBt1Vbuwx+M7WBkpfU0oaPFtS6ic6fHpMUI3GRzXQwhdB8MwO3bs6Jdr32QyFRUVRUcC69evH86lBt0+omSgkxTDMOvXr+dHb9LS0rZv3z5Os3GCnL+dqWl1Ho9+otVIstLpwnTFnLFfPzoC0ckG/wpMwOJlHU6/xcs6vWxvj681em+/JgD0QP8tvH4DHfSVSVYTiQlYgpyfX3YS5HzOgAUG+wq8sM1FWJyExUFYnESPm3B4CMZLODwk4yZdTsLplHpBOtiJnAL8d69g77wFFmRGb+eDBymlllJKqVAt7WveXC/RfWzIZHFf7HKxXN9rUFpMFaTS2RrZGNZjQ6D2fOD40UjFOuH8heJ5C0Aw2glUCN30cKQCoUkpunZvRH5+fvTC00GP4SeaR5qR7Ph8wv7ICENpaWkkweVQwwsAsH///sgd+XQ6aWlpfP1jPrFpQUHBwOlPkdCivr6en3+QTDYAACAASURBVAplNBpjpwSdAiLlI4qLi7ds2TLSJcLR66pjaGdqzK6zkeaIShGPK/6NvtNv9bG9zgFfZDjjGzw+zOi30cc6B44JROs3NqIU6YWkBIYYTOAFOb9riB942OYiulxENwNmB9npIK29pLWXtFgom2vQ42Mg5ITojoeSv/govfz2yEbhYN/xplVvddW0Xx0LUkuFM1KVGerRr8YGgFBvr79qH2ez8k1Sp5cU34nznRAaIwwqEJogDMM0NDTwn+vr6/kHa5qmi4uLI+sKousl8wkx+e00TUePCfCP7wNvET1Bf9DMkrxIKhUAWLly5aCjB9GXMplMA+cmGQyGoqKi8Zi9s2PHjsrKytLS0lWrVsX94hOJL2sdXQChsrKSYZgYPzT+yZvl/IzfAgDOgDXIea/7xDwoCaXK0yzXyYwT/2p/7Pifg93bMuhAx8QIe4OE2QFWJ9HpIC0OwuqkGi+TXR2C3tD1T76CUmsFqlSBKkmoSRKoNJRaQ2m1lEpNqdWUSiPUalW33jZ+X2EKaOpyn2rr++NPVoiNOtqgGtvfM8f5T54Inj0FHAcAIBKLFt8umjFzzD1FCOH0J4TiwWw2NzQ08O/d+YDBZDLRNP3mm29GHt8jU1/6OXHiRGTZbmVlZeRzP9XV1dfNERS96HmoBdAKhSL6HXl+fj4/oGE0Gvk6awaDQaFQRD/4Go3GoUYqxgOf6WXLli0VFRWlpaUTkxkpjkwm0+7duysqKvhozWw2R8KzoUZjvKyjseeQ1VMfl7XUBkVhOj1nUk+s52f4RH+FfmFGv+OjlxYMxHL+fmMjPb5WON4C3QzlYiVuATjcYcYNvS5we4krleCI0zXD6SpBApWcJTJMExmyxOkZovQMsSFDnJEhTs8Qp6TAWGbn3PTaerx8RCEUkEtzk5LlojFeMNTZ6T+4n3P0/f0IcvIkS5cR0kHnoyGERgyDCoRG7OTJk9ETjQbWbOIxDNPR0XHd7EbRD81FRUXRtZAjtY0LCgqiD3v11VcHxif8pKPoS+3fv3/gYf2e0ScyWhimbdu2bdmyxeVyVVdXr1+/fvXq1U899dSNnySKYZidO3fu3r07evEJn1A/xll2X2uL48RQaV757D38fHqIWokbw1SdcD8wzOjH3WDyXbrkbbnkb2n2XWxgHT0hxsF5XZynl/NYOd81Bw+ceXbdZ39KkxxmfeLsQsk0oyQ3T5I9TTptujg9Q5IRazk1GjUL4zvWYgcAAUmU5CerpGMqYh0OBAJHPw2aavkmIVeIl66gsq7/PxRCaPhw+hNCw3Ly5MkTJ07wCY4AwGAw7Nq1i9/VL6jgI4GFCxempaX1m+jCnwsDAgDUD/+AHh3wRCdKGnThOAAYDIYNGzZEzzXauXNnZC9N00ajceHChZFQbfjMZvPJkyf5O0ZPS9u4cWNkxcjGjRurqqoipxQXF/PFuYZaO9Huqml2HI+eta+RZOlkBXxGpkk9zhAXbE8P63aHPO6Q28253SG3m3W7OY8r5HaH3G7ptFxfyyXfpSZfy0V/c13Q0jKWe5EiECjTBWqdQJ0kVGsFag2lUguTkkVpBnF6hjg9XTY9N17fCw2H3ROoaugKcWGCgBV5ySmKMU15Yi+3+Sv3hX1evimcNUe8+DYsQIFQ3GFQgVAsJpOprKxsYG4io9EYXY/JbDZ3dHRMuok6Nziz2VxWVsavI1coFJHfAj+CMegp0YXkSktLB12qXlxcHD3HzGw20zTNP/rz617MZnNaWlrkOgzDrFy5ctDbDVzObjQaV61aVVJSMlTQGOT8rY7j7a6a6Dk8BkVhnmbZ1BthCHR1BW1Wv9UatFmDNmvQ3h3oMLM2C+tyhoPBcDAQDvrDQX844At0NsfrpqRMIEqbLdKnkzIFKZOTcrlAJhfQNCkSk3I5pVCQclogl1M0TanVlFoj1uOS6BuL0xfcb+oKhjgC4NZpSaNekx0OBIJ1tcHzNWF33wJ6Uq0Rr7hDoBtyphxCaCwwUkeTldls7uzs5D+YzeZVq1ZFHuNMJtOmTZsGvsnuFwls3749+k02j6bpDRs2RF6K81Vdo69QdEX0WfxShDh9s0kpUqbA6beGOB8ARJJvDoeQlGokWUqxLvoNvcFg2LJly+rVq8vLy6MnEfEP7oOOVERnmuo3J62joyMyrywiEp8YDIZ+F4wsZ6dpWqFQRAr2wZXBKJqmo3sVHWAMauDCCYoUZysXZasW32hrqYN2e8jjCbn7hgU4fpTA6w65XCGPJ8wGY5wbDgR6973PdrUFu/r/tOOIFIEwbZY4PUdoyBSnZ4jTM0SGDEl6uiQrWyCTjd990XjzBEJVDV3BEAcARdma0UUUod5etuZMsNHUVx4bgJDJhTNni+aPeJQSITR8GFSgSaasrGznzp0DExbV19dHXj/X19cPfOIEAJPJZDKZIs+dFRUVA49hGKaysjISVJSUlJw4cYJf0lBUVDTSNKDjamCO0bhk8R80dWmkpFd0qDC6rERDiawrUIj0SZKsSA3ggSEcTdPDqULNzz6K3sIvoI++WuTvpN8fTFpaWvTvuqKioqOjYxRlJZiApdvb1uNrZQKW6J8Vn5opLlUjRifk8fhamv2dnb62lkCH2d/WGuy47G9vDppruVghw1gJk9Oo5ExKrSWEQkIoJoVCQigiRCJxdk70YQQlFMhkAoWClCmovg8ygVwukCsEMplQqx3HLqLE8bNcVWOXn+UAYFaaMkc74viQbW0JnqsJmdsiW8gUneiWOdT0PCDJePYVITQATn9CNxx+9UJHR4fZbOZzKL3yyiuRcYCSkpLod8YR0WWbGYapqKgYeFi/RQ6RBRL9RA96JBb/+p+vugVXEozC0IW3biiD1hMYipd1+KLKEve7SJIkG65UWe53wNjjKJPJVFlZWV9fbzAYaJrmQ8exlKWz+1p7vK3dvpZBay9oJFl5muUTv2TCebLacfig8+hhf0tjoKOG84wgL+qgSJmAFNKEUEpQYkIkIYRiQigmhCJCKKQUSipFL0zRiXR6YYpOpE8VJScLU3Si5OS4fBc0JQVD4coGm8MbBIDpyfIFmerhnxsOBoP1dWxtDee88s8ISQpyckWzC3GyE0ITBoMKdAPZunUrn8i/3/boRbqVlZWRmsT8o/+kXvQcx7oEEyZSRyw6j6dSrKdICYwh+5CXdTABi9NvsftaxxI1Ka4NM4SkVCnSAYBEqOKHPkZ95UFF13Tr9rUOWi6NIsVaSbZGkq2X50/kwgln9QnHp4edhyqZz97tl/6oH0qTLEzLExuyROmZQp1eIJcL5HJSpqD4D3K5QCYTyBWUXE5pNBPVfXQT4cLhqoaubncAANLV0tunDXcwKhwIBE4cC5rqIDIrTywRzrhFdMtsAifCITSxMKhAN4p+q2/5ws98hbUYKXQmC/7Rkx9z4COH4ZcQjpBceV5XivQUKYYBOUbt3jHlwIkYmLp04t+sMwGLM2C1e1udActQhY1HgSLFWkmORpKVJM0caYARiR9CnK/b1xKj4jIASCiVVpKtlWZpJVkTGUg4jh9zfnrYeeQA8+m7XKD/XlImoBetFqZnitMzxYZ0SUamyGDA1EYoscIARy52dzh8AKCjJctzk4ZZ3iNw4Xyw+ngkrROpSRLOmSs0FoxfVxFCMWBQgSaUyWTiM3JWV1dXVlbSNL1r165IwLB9+3aGYSZvFME/dMKVFQge1uFle0c07MC/ZY+MAPAP98KRzCOawiI/3mhOv5XlvNFbmIA1GPVaPsajf2RulUaaFYma+ClnfOQQWUAyzAhQI8lSivRaaZZGkj2ua6/Z3l5f+2V/e7u/o93X1Ohvv8zabWGOY7utvqbT/Q4WKGT0rQ/QS5arlixTzps/fr1CN6ceT4ALA0kASRD8f+XikS3XPNHS22x3A4BGJirJTxaQ1w8pQhaL/3AVZ+/mm2RyivjWpYK0tFH0HyEULxhUoPjjI4fo8nBms/npp58edPF0eXn5WOavJ0SQ8/f4Whm/pdvXAjEfWwfFDzjwFc2kQjX/GhvrEoy3K2umW5wBy6DrN0ZEI8niYz/+NziuJed8bW3M6VOuMyfdZ08GLl8KdJyLPZcJMJBA48wXDFmZgNXlt7n8bj8bl2sqJcISY7JIcJ3l1GGv1/fZ0VBjHd8kZHLRwltxdAKhGwFmf0JjZTabd+/eza+r5v/Lby8qKnr11Vf5z9HbeXwdgNWrV0+iiMLLOizuBqunfpgz/vkX4fyjJ0VKlWLdVK12PCnQIj0t0ueoFgKAl3XYfa12b6vd1xI7wOAXkERHgHFJsRVbwGZznjjuOnvKfeqE++whtqcrxsGkBIS6GUJ9plCXJkpNE2dmqW5fShfOHdceopuQNxiyMgFbXAOJCJlIsDwv6ToRBccFztUETp2AYAAAgBIK58wVz52PZewQukHgSMXUF6n/FV1I2GQy7dy5s9+DPk3TfN2uyJadO3cOrPtmMBieeuqpyNrofoWEI/oVhSgvL3e5XPyi6jFm15lgVk+D3dti9TQMnMUUWRPMjzNEhh0m4LkTxUuQ8/f4Wpz+vrEmfsrZxId/vUcOOz497D51wn3mSLCrfdBjRGnTpTMXiNMzxRmZfcWeDenSnJxBD0YoLtp6vFbGb3P5XQMCCVpC6WmJ8HpjC8MxLUkmEwliHBDq6PAfquQcff8IC3KNksW3EXL52G+NEIoXDCqmoLKysurq6oGDA9FJVwet+8Y7ceIE/8FsNq9Zs2bQY6LrfJWXl+/YsSOyrtpgMPCRw41QXnrQkgsAEMm2BAC0uG/Rcz9e1mF1N9h9zf0m00solV5m1EqzdLJJExehG1PvkcO9hw+4jh5yHtk76AHC5HT5vKXyeUX0/CLlgiJKPYIkmwiNBcuFW+yeRpuL8V0TSwgFpF4p1tOSVKVYKowVBsRFOBAIXW4NNjWGWi7xW0h9mvj2ZQJMT4zQjQcHDScfhmEaGhr4MtImk8lsNhsMhm3btvF7zWZzWVnZoCdG510tKiqqr68feMzChQujjy8uLh6Y4BUAokcz1q1b1690caLwk+aZgKXH1xr3lKwaSZZOVjDBKUHRFBMOBp0njvccPsgcqXJXvz8wOxOl1ckLl8rnFdFzF9BFC0UpKYnoJrqpuf1sY5e7udvDl7UGAJIgtHKRnhbrlWKNTDS8zExjEurtDbW2hNpaQh1XR+0IBS1edBuVmzf+90cIjQaOVCSG2WymaTo6wVF1dXVVVdXAB/3S0tLIK//q6uqtW7cOutx5165dkZhh+/bt9fX1NE0XFBQAgNFopGl6UhdziIEvNOYMWAcOKYwdRYp1sgJ+UAKnM6HRYZ1Ox2dHmepjzGeHmaMfDDxAkjdfeVuxomiRctFiGT4wocTpdPqautx8alceLaGMOjpTI6WGkZFp7NjLbaG21lBby9UadleIFt4qmrdg0LMQQjeISTNSYTabT5482e+xmJ9yM4qr8S/7+22MzlYEAIM+4peUlETf0Ww2V1VVDXyXH12S2Ww279ixg2EYhmH4qm0R+/fvj9xx06ZNg44JlJeXR4KK+vr6fhGF0WhUKBQlJSXRP5mNGzcO/rVvYNF5PPlMrJFdUkotGzA4wOf6HHTBNF9ojE/Jyi+PHnhMJEnrUPOjABO5ojHwNDY4jx11HjvqOnHYd/EMXPvqRiAXyufdo1i8RLnoNtWixQKFIkHdRAgAgOXCzd2eJpuLubJqgiDAoJLmpchTFBPxMoVtagg2NYXMl68WsOO7oaApQ4YgI1OQkUmIRBPQE4TQWNxwQcXJk/+/vXv7jqu68wT+23ufS91vUkkyFnYDsRFGpLkEkhaJIUlPJosmk+WOWSt0v01PRNaa1/4D+Af6HWomvNHOamgYA/FKejrdbrkRGWhBMpFl2QYGG9lIVZLqfjm3vefhVJVKF+tWutn+fpaXV+nUqVOnbNU5+7f377f3x35iT6VS+dnPftbZ5l7RIvc9//zzL7/8sv/4ypUrL7300uqm+fHjx1999dX2oV5++eX33ntv9aGeeeaZv/u7v/Mfnz9//m//9m9X75PJZDqnQL3VWV2+fLl9KH9BhjU/bLlcbp/V888/f+bMmc5aal9nZtGPfvSjSqVCRE888cTtUutctuecjgEEP2yg1hoOm5mMNU8bT7XUXmhsq0sm69zEXK4b8mo1r1p1KxVZq3rVqlevy3LZrVa8ek1Wq161IqtVz39Qq8pq2WvUZbno1cqyVpC1rD//KTeI6QFmRJgWYLrZ/GOYXDeYGeCaxnSD6QbTNGaYzDCYrnPDZJrGDIMHAiIU1pJJLZbQk0ktkdASST2ZFAepTLP+xRe5//WP5Q8/qPzHb7xyZcWzWiIV/dYPY99+Jv7NP4sMP7IvZwjQyfbkbNFarNlfLFRd2Qx8DY3f1xP+Wjq8K/USrusVCyqf94oFmS+oUr690MQSTReH7hGH7xWHBwWWbwe4rRyIoOK999579913VxcWP/PMM+1O+uPHj6/ZfO8MIW7evLlmZ/+VK1euXLnSPtSax6HlJQe3yhSKRCKRjm7FQ4cOrXm0zhrlJ5544vnnn/c/WrtiIRKJfOMb31gxvLDhCEM0Gm2XRx8E/vCCv0o0EfkLRW910YbVtOVDBLdadyxi9PcEjuzBQmN3jMbMjJ3LOtmsnZuzs1k3l7Wzc+5iTjbqyrGlbZFrS8dWdkO5DWVXlFvbcD2EzZM2kd2g6s4dkYiItFSfiKW1ZFpPpbV0n97Xb6T79HS/2d+vp9NGuk+LxXbkjeyFhcL7/+7k5uzsnJOdc3NZZzHnzH3p5j6Tztov0RKp2Ld/FPv2yfg3/yw89NCOnAZAl2q2d7PYuFGo5yrLrqvxoH4sHbk3GdzM2nNtqlyWlZW3XXdFjq7jyMKiLOTVqj3beLpPHD6iHR7EAnYAt6/9r6m4cuXKmmW+jz/++Msvv9zZ7PaXQVixm18w0P7RD05W7BOJRDrf4ubNm+fPn/frDdru1JKD9iLEnfMd3Yq/lNuGNrm68PrakYM/UBAz+zUewDIOG/IqFSs752Szdi5n5+achXla9ytsffH/Kp984OY+83a6Nd+J68TDKR5M8FCMh6M8ENzqEZRtKddRtqUcSzq2cizl1JVTU06l+8DGGPgT876Hot/69vq7eaWSW8h7hbxTzHv5ea845xW/ulXAcCtasjf6rf8ce/pk4unvIJCAg6NiuTOFxo1CPV9bNj/AVjOdVL3u5bJeNivns14uS9a2vp+GKXp6eaqHpXpEqkdgQgKAO8L+BxXUyvA5dOjQgw8+6E9Lelsk9uyXxY5CgvYoga8dFexIu3+rViza4Gsv3eDDAg7ra8YMuZydzdq5OWc+Z89+5eayTm7Wyd1w5z9bPV/QtvEA6ekHtb7Dwgw0U490vZl6pGncMJmuc8NgpikCQR4Oi3CYhyIiHNbCYR4KiXBYhCMiHNY6ovrd45ZKbqHg5BedQsEtFNxi3s3nvVLRWZj3Cnm3VHCyX7kL193C4h6cjIjF9d4jet9hrafP6B/Q031GX58xcI95+HD4waE9OAGATcrX7BuFxo1ifcXMsBpn/bHAYCI4EDPXX2hCWZaXnZPzOS87JxfmVa265ZPgnMcTLNkrUine0yNSPVhfAuCOdCCCiruEt7hIjsOTybycXXOHfH0pWvALkVuPu00o2iRtc6XJ/urC1BpeQEHzJlmzs24h7xYKbqHgFBbdQsEtFuy52WYizez1bccMneUKTm5GREI8mOKhGA/HRCSmJVJG34CWTut9/WZvn97fb6T77uAV06y5OXtu1llctGdvOtmsnZ1zcll3fs6r15RjK9tSrq0cW9qN5mCIW5V1RUTcJBEbFIk+LdGjxZMimdLiCS2Z0hIJLZHSk0ktmTT7+s07cUgT7iSLVXuh6sxXrYbjLVSXXVN0wQ8ngvfEzXvi6w0nevPz3tysNzur5udkubTmPiwc4b1p0dvH02mmG8Q5CUGcM875DuUcAsDt5YAGFevMybN5K3rxd8OKeYrauKPMihLlBi9Vg1UeqopAg7PWv7RtqFrIq4VlLew1QrIW9pyA3JHzSbaGCNrt/vb0R2j676rGtWvNrKT5rJ3LuXOzdm7OWch5+ZxXzG4ji6bN7xTX0vfovf1LneKHDqOBCwA+25O5irVYdRaq9nxljTHqgC4GE8FD8UB/dO2BYmVZ3uxXXjYrs7Nedo68lYtnExELhnhvmqf7RDot0v0sENjhjwEAt7kDEVS8f+O1Xe2JNyyuW5x7TEgiybgkrogk40RMEpOMK2KSGJHkJBkprhQjyYi4kowUJ8mU4qQ2ql7rndNDNR6qCsPe2pTenlC1kKyH5WKPY5uKiDSmh/QEESnOPY3Fw4ek4ErjJLg/PuC/8ICECqpSkeWSVyqqclkWS9SoEefEOXGNOCfBORck/C3C39L1ezIWDPBgiAVDLBjkoRBpm5p1QNm2sixlWeTYzZPhnIRgnDPB21uYrne+qv7FF36Js5Wbc3JZ189NWsi5c18681f8fu7tEbG4nv4TPX1I7+3X+/r1vn4jnUbMALvLdZWUJCV5npKSpKekJK+5hYj8LmcSnHHBWt8RvxN6k1802G22J2fyjYWqtVC1K9YaMYCp8WTISIWNgaiZCi+bj1XVarJa9edvda5cltk5WcyvPgLvHxD9h0jTRU9K9KZZGHMfA8B69j+oKNtz4zde295rNYdxlwmPCZdpDjccpltMt7lhM91mwmOmzQ1rL5bsWYcbNt2o4UZNpXEiIkXhukalKitV2Vq9QRtigSAzTQoEmG6QGeCmSWaAmQY3TTLMPZjMW9m29OOHUlGVS7K4w2tXb5Om82DIY+S5rnQ9z3VlrSbrNdWoy0ZdNRqqYclGTUlPea7yPOV5ntVQti0ty3MsZTc8u66curRqnlNWjtIiKWLcrc5vIyVJhANa7/1632EtnhSJpNb808qiicf9SVGNnp7Vr1WWRa4jHZcch1xXSbfZ2pMeSSk9SdJrb1FSMdPwgysKBnkwyMztlKwoxyHZOqYnld+4bG0hoZGuc00jXWe6TmIX5pqEnabqdVmvU6OhalVZr6l6XdXrqlLuXKJ421gozEyTmSbpZusqZDLT4GaATIMZATJNbposuOWSfVif48mZQmOmUJ8rrayQ1kUzikgG9WRIDypXViqqVpXlsqpWZKVCtYrfAbTO8cXAId5/SNwzqPX3I4AEgC3Z/6CCiG5U/lh3VjZMe1Q/NWxq2KzeoLrNLItsh2yH2S5ZNnMcVq3vy9muR2g8HmfxJE8kRTLB4kmRSKzTAlPlslfIy0JBFfKymJeFgmocvA+1T2TDkrbl2raybc9qKNuRtuVZlnRsz7al3fBqNa9e9epVr17x6kWvsejt4vxGS7hOIhgTkbSI92qJpB5P6omUHotr8ZgWjevxmJFMkd/vqzbMalPKdclxlOMq19nmPCoraDoPBikYZIEgD4YoFCLXJdeVjkOOTa5LnqschxxHOY5yXXK2VcZhmEzXWTvMYJyag1GCeMewj+DEOQuFxD2HRTSGNsqGlOOQZUnLIstSdkM2LLItZVntjRv+UinbUvX6dqppd4lu+OEHM0wyAjxgkmEy0/SjDjID/rN8T8r9b1+uVF8VG9fztc7lrokoFRAD3EmRHXEbeqNO1YqsVFStIiuVFQvJ3ZLQRDrNDx0WA/doAwPoMgCAbTsQQYV7/Zp344as16heU426qte317ZmwZB/r2KGSWaArAYPhSkUYsEQD4WYGegY0+d+02flmH47K6CVGKCkR+3EAKWW0gBE6wit5JmdajApy2o1AR1yXeXYK3/saHOQZUurQfZeT/RERKTpPBplkRiPxVk0ymMxFouLRGLNfb1Gw6tW/YXSvFrNrVTcQsH1i5Xzi24h7xbyXjHv5he8Ys4tXusmoWjzuEHcjPBAggcjIhhhZpAHgiIQEIEQDwR5MMiFEJwLzoRuiGBICwX3ZqajO5Zh8nCYhSM8HGWRCAuHeSTKwmEev0vnEfYWF1Vh0csXZKGgiouyWFwzl717LBhikQgLhXkkSkbX06+5bvMqZDeUbS3lE3bJDDA/RtV1tmHTlonmCEnAv+Ab3Aww0x82MZgZuDMiWE+q2a/m53P5cqEsbMt0bdNumK4V9OyQtHTLInez/+wsFGaRCAtFedT/TQizcJRHIiwU2tWPAAB3j/0PKmSpVPuH19fbwzB5OMJCoWYvVzNmMMmyeCjEQmEKhXgotAdpPweZajSUbSurQe4ONEqk1ZCOIx1Hua5yHOU/9h+4rnI9KaV0HX9BZbdSlrWarFa9WlVWyrJe8+o1r1KStbKs5mV93qtut0h5XSIaEfF7tESvFk+KRE8zv8ifqyeW0BIJbppM17muM8NguiEMgxkG1zRmGFtdiVk1GrJUUqWiVyqpcklVSrJYbHYGa3ozquSCNetG2sHqJvr8hEa6zjS/LaUxTSddJ03jftNK02lz6XtKKqpVZL2uGpasVqlRU82Ml1rzJDln7bPyy1qa8bAgIUg3mKYzXSNdpw1PWylybHJc6Tjk2uS6ynW6HG9hoTCLRlk0zmMxEY2xWJxHo3dUc8d1vfl5r5hXhYIs5FWhIMtbmEaCBYLMDFDA3PiXStN4JMrCEb/VyCJhHt2jqXhUvb6iv0NZNlkNadlk1ZVt+6MuezaKstTNpJtkmiKZJE0nXfe/XKJ/YOfvGv4F03WUn74ovTX3cjzleNKVynM92Wgoy5JWg1kW2RZ3bN4aYVBKBRe3Vm3IAkEWjbJgiPv/9eEIi0R5JMIiKIcAgF23/0GFsu3pnzxnz3/FjQAzA8wMskCQm0F/eIGHI8xsJVqsjzFumNww/L4uruuBIysnzfQsS9m2dJzmSluOI21bOY6ybc+2ichvdzJNY4bBdd2fvL/VKtW5ENK2/SNI21aOpVzX39I81O70Mq6gLEs239HyG/qq+cBWji29tW9jgq6OWAAAFkFJREFUy47QqMuOuTXtG1f34LTXoaX6tHifSPRqiZSWSGmJpEgkRDjCQyERCmvhCI9GRTAowmEeCmuRiAiFthoVwB5bURlCpJTnyVJJlUuyWFSVkiyVNk7P0HQejbFolCd7tMHB1lZGhtHMolleT3/QeNmsN5+VuZycz8r8egto8ESSRWM8GmfhcLMawQy0EoQO+sfcDtdVliVtixqWshvKsmSpRJYtHYccixxn4yyv9mjJjiQNrkM3mgE546paJjPANJ0ZOmk6E1zZdiuQ2G4m4VYZJg8GWTBIgTAPBlk4xEJhHo2y0N074gcAB8T+BxXVS1N/+N7D+3sO0D0Ri/NgQoRjPBQVkRgPhvyF0kQ4wsNhEQrzcFiLRHgwLMJhEYlo8bieTGmJhLiTeqNhK1StJstlWSzIckmWyqpSVKVSc2hl01q90QHGN5pSTNf9IhPmV7QHwywYZMHgjnVXS+ktLnrzWZXLefM5mV8guaplbAZ4LMYicR6LslhMxOI8FsOkOl1qT+mmLIusRnO0xLLIbkYdzRytSnm/z3Q9HheOMFxNczTT0wypG56mSzMQjkV6e2N6OMJDIQw4AMBBtv9BBRFd/u+jjS8+VbalXFc5lnIs5a9L5VrSKcvaxl3vdxseEkyLcD3A9CAzAkw3mW4yTWOGyfQNWkiMMabrTDebQzGGYR5dOaTDW8sqM8MQhsE0nRmmnzjEDEOY5tKCyqGQhu4x2EGu65VKqlyUpZIqlWSpqCrl3Z1hTGgsGGThMBHzp9jqmGtVkpTK87rphOY9vaJvQPQP8P7+PctE2neeVFIpqUgqJZVS/gNJnlJyh246nlSeUkox/y2kbL6X/6ZKkaeU6thetT3NsYXncv+P9Jj0mJRcSZKSK8WkZEoyJf2NNcvlUnIlmZJceVyqhhHUPEfzXCFd3XNJKU8Il2ueEB7XXS6k/yMXrhAe1yXnuuAaZ7pgGmea4JpgOue6/5gzTRMiEBABU4uENc7WX9kaAOCAOxBBxU7xGg1l2+2MpsaX1zd+zZZw7udWNbPz9Vaqld/UxuSJALtG2TbZtrQaZNnKqkvLpmbKfrNimDa6lKlGQ9VrezHBmmGKvn7eNyD6B0Rf334lLylFVdutWG7V9qqWV7Hcqu1WLdeVd841f+/pgmt+hMC5JpghuKlxXXC99cDQmCG4IbguECQAwN3ljgoqAAA2pPwFTGrNQnZVb8h6jep1ItVcn7FV0c6aFe2t2XI3wkIh0T8gUqmtnpLjKcv1Gq7cxgU5W142imK5Xtlyq5ZXs/eixOuA0zhjjHHOOKOwsbIwjzMmOHHGOGOMkWg94IxxzgQjQ3BD47pghsYNwSPmnTCjFADALtn/oMKV6n9fylZtl4g07l/9l67v3L++M2rfGHhre2jVHWJ9rHWo9hFa76V4x3uJpZsKccYYY5yRYEzwfV5EDwBuU+2YoeFIy/UsVzYcz3/ccGTD9bzdHz0Im1rY0KKmZmi73n2+zsVWtC6q3fMv15wTo+bleukyjss1AMCe2/9+l4bjVVs9aq5UdJsPzeuCC+4PjjNtrccr77X+LZar1dvbN2PW2qLhTnkbcjylSEnZTvVWUvmZ5aSUcqVyPeVJ5fp//IkmpXKkcj3pSSU40wTTONc4E63fJf+XSnCuiU1OObseVyrXI09K/xw8qRxP+ufgSel6G3c8LP3S8pW/w4yRf9rNHl/B/d5fQ/A7rOVnezJbtv04wf/b9mTd9vYmZmjTOIuYWtjUwqYWMbSwKSKmCOkau6P+sQEA4MDZ/5EKIrq2WC9brlxeVNeu8JOKpGw2xdrbPaXq9t1bwN0c019qwy215JIhY/1E3tWdiGFzy0uorjgC62hNtiMfx1PLajTb/4M7VKcplfIkazfTd/wtlCJPKldK1/Nb/NKTyvH8Fr/0pNrLluKdKqALQ/CwqSVDywoPGJHGueBMF4xzCuoiFdr1hWi8Zkwl/Uiv9Z9OrR+l4ynbk7Yr23/70dc23svQeEATpsYDevPvgMYD+saBFqNmck57FHerY7YAAAC74UAEFXumc0KSpZZox4Qk7elElFJeK8JRzX7lTf1TdXbxOs0e6GYbxfE2mHkd4IBYGnATTNtwqlaiteK65rdsz8qCRXskxx/VEUuP/U/hKeW1xoL8AZnOxzv79WzGDDoPaCKgc1MTAY0HDW4IHtBFUN9yGA8AAHDA3V1BxUGwdmCjSErlF2m2dmgGNqr51FJg0zmAI5cP4HRM3bj2dnSub0krxYh1pBu18tn4prJ32tnkKwd2GONMdbaDBWfm5pLdbU+6nvKkdHZotESwZrNbND/mzifK2J50/A5+T9mucrxmZ3/zgecPBx308Nuf9NNoTvLDm5P8aNzw07o0hpgBAADuWggq7kbLAxuqbn2WmM7IR64MY5QiYoyJZsXIsvJ3xpjgiu9EfndnHX+rNGVlChbcvmxPeq28I7f9GyvJr0VZmSGpWnUp6w5EGBoXnOntAhWxbCijXQHlF6u0o0dNMMGZgelBAQAAbg1BBQAAAAAAdAV9bwAAAAAA0BUEFQAAAAAA0BUEFQAAAAAA0BUEFQAAAAAA0BUEFQAAAAAA0BUEFQAAAAAA0BVtB4918+bNHTwaAAAAAADsl3vuuWfzO2Odij2VyWRGR0f3+ywAYFfgCw5wN8MVAO5ySH8CAAAAAICuIKgAAAAAAICuIKgAAAAAAICuIKgAAAAAAICuIKgAAAAAAICuYPYnAAAAAADoCkYqAAAAAACgKwgqAAAAAACgKwgqAAAAAACgKwgqAAAAAACgKwgqAAAAAACgKwgqAAAAAACgKwgqAAAAAACgK9p+n8Dt4eK5V96f4cnh0y+MpFrbcmNn3po7+kLHlh0wO/76O5NV//HytyOiqbOZC3PEiIhFh3/yYudzS08RHf7+6F88sPywMnfhl29f0lceEAC2z/9aVTq2rPpiAsBt4+K5V8aLX9+Nr3DntWLVVSI3duat6fLat2+/7eE/PnJy9IdDmzwmwP5AULEF+cl/+vDYT59K79bxL5575f2Ze78/+tcPEBFNnc288YvCd//muePU/PGCPfzC6EjKvwb94xnyLyKtgMF/imbHX38n88uFU+3zXLpgJXfrxAHuXh13+tzYmbfezFx/9NQuXiUA4PbiTp99bWzuyMnRvxoiIrp47pU3M/lW8JAbO/PW1fh3R1887j/128yvqOOpy/T106PtG33mjcVmz+C6xwTYN0h/2qzk8DND0eIfP7qyS8eXuQt/mOFHTrYvCie+ORzxZq5+RkREs+OfZKOPfL/ZE5H+9p+fCJcnP5wmIrr00cVqdLj1FA2MPNZPxWtXF/0f3enx2mMvjY6eGopi6XSAXZU++eJfDkWLf/jn8cX9PhUAOBhy45/MisFn24MMDz35cIRuXJ0mInKnx6fLg882uw7p4edO9tPMx+OLrafif/rnzVs7T3/nTwdl4dr04kbHBNhHGKnYvPTIYwPTY+d/PX28cwiyZY0RzNbIQztOWD9j6qEfjn6n84neVISovJCjB9K5K9cqPN7bfpanhwajU1c/v0JDyf4nf/6zZd2i6WSUrhbmiVJEpA39+IfdfWwA2DT/KjH54fSIf5XwOxT951bnT94q5wEA9tfy5MbDnQMLc0dfOFH4Bz8raTN5Ryf+/KWTHfdonu4JkyovLhKlLn/+FYt+vSOJ4MTXBsfGr00vjozInqdHR3s7j9OTiKqZfJ4ote4xu/rYAN3BSMUWaEM/fnpQXh/71Wcrn2mNYI6Ojo6O/pfhwm8zv/qs1evQ7jyQuemZ8uDjt7j+8HTviicuf/4VUbQnTUS5fJliic7rS7onrmRxfpHSvcsTLWRueqbMlu8MAHuE96QiROVFv6/x7Gtj2vf968Loqf5rb7zRHMNYdsV4evDLfzmDwQ2Ag2Pq3ben9OHTa35D85Nvfnn/z0dHR//bqYfaKQO3tvIe7U5fnSMWTaWIcgtF1tldSH7kUM7niXrTK27iuSvXKiyaTG5wTID9hKBia/xBxo+XNwBWjGAOjDzWGsE88bVBdePzZsbUpY8u6sN/tsn+SD8bKtmx/yavF5c+ulihw7cKXQBgr+TGP5ntSGhMjzw2kJ/84LNVV4yHnnz4xNGe+dy+nSgAdHKn2fDoS+1xxQfvP+Q39H3J4dP+OCRPDw1G/fGBzcuNfzLLosNPtfIdNtkD6E6PT5fZvY+teWtfeUyA/YL0p63h6e98b/j6O5P/9OGx/9TeOL9YYdH7OkYwlxKQHrz/0PtjVz+j4w/Q1Kcz8aNPbrKtP/Xu21PV6CM/6biAbGZkc3b89feXFWYAwH7J5cs0N5bJjHVujC3kKLz8isHT33kahd0AB4Y29NADqzKgduTIF8/943Q58eippVt7qZWrvA6Zu/APY3OdRRTrHxNgvyCo2LKBkR8MXXvrjx9dPtbaslAo32pnbehY/9iFj8cXj6auZqNHv7OJpoN/IatGH+nI1FxWJkFEaw6b+nPPrZh1DgD2klxYrBAdSaVk7mKV2Jrfx4sf3fKKAQAHgD9Le/zRU6NPpf08xu6P6ZdRJTpmh0v3xNXl4vwiHW/fyhcK5eV9lM26LDH4bGsqyPWPCbCfkP60DemRxwa8mclWkeXKsdHlJRAnvjkcKVyb/uDzrxJHhzbsSHCnz/7Pty/VB5/92bLar/TxoxFZnG8PsvqFE4fvb19icmNnXn1/JvHoKUQUAPtoKQ9hndSI1VeMC+feR00FwL5a6qq7eG4sG33k9OgOttSnzmbemi4Pfn/5MVddCqY+nWGdTYXZ8ddfG5tLDp9eK6JY+5gA+whBxXb4FdsdP44MRWfOn2vWTsyOfzJHSwXZfceOhMuTl2bu9bfI3IW/z2R+cW6NqWn9y8eRk6OrLx8DI4/1lf/421aV57//81R1KYFy6mzmrcu0s1dAANgqv9ewPQtk+vjRSH7yzV+36jgvnnslk/nlh7nVV4x/ujTz5ae59S4OALCrOrvq2tXSROT3FGziALmxM69mMqvncfG/1/8+P/jd0VWTvK24FFw8N9bZeLh47pV3JsWjp0ZXzxi5zjEB9hHSn7bp4edOfpq5YDd/Sp988S/pzFuZzHkiWjFBZHP61/ix9b/57vRZfy3t68uTsFvpEyd+PEpnM29kJv0VtdvJUbmxMxfmiFF58s3MZPtV7XnuOrJCGU2+mZmkW4yiAsB2dH5hxWBzESvfwMhf/9fU2ddaO4jB742O+s8uu2Kw6HBrfau9PnmAu1xn4US7AttPcv5tJvNbIqLYIyeHI2PXF3L0wHa67abefXuqQoxmzrdaCETUvhGvaDwszUHvV0gSlX7/dub3S0eLPXrqp0+l1z8mwL5hSmFNNAAAAAAA2D6kPwEAAAAAQFcQVAAAAAAAQFcQVAAAAAAAQFcQVAAAAAAAQFcQVAAAAAAAQFcQVAAAAAAAQFcQVAAAAAAAQFcQVAAAAAAAQFcQVAAAAAAAQFcQVAAAAAAAQFcQVAAAAAAAQFcQVAAAAAAAQFcQVAAAAAAAQFcQVAAAAAAAQFcQVAAAAAAAQFcQVAAAAAAAQFcQVAAAAAAAQFcQVAAAAAAAQFfu5qBi6mzm1V9P7/dZAAAAAADc5nYlqLh47pVMJvPG+GLHttzYmVeXb+nK7Pjrvzh3Zc2n3Omz/+NM+512+H19Mnfh7zNnPsxtsNs6JwkAAAAAcMfQduGYU5/OxKPRUuHa9OLISGoX3mCHnPjx6In9PgcAAAAAgNvezo9UuNNX5yj22GOHVPn6pxv15W/bQqF8q6fmFyu79a5btM5JAgAAAADcMXZ+pOLy51+Jwe8NDbmXxy5cu7r4VLpzrCI3duaN6TIjouTw6Rdawxgyd+GXb1/yQwEx+OzfPHfc337x3Cu/o+91/jhe/PpPXjw2eeat6TInOp/JnO88jr/P+zOcaPLNzKQYfPZvnksu304sOvyTF/0XTJ3NXDBPvvTDoWU7EMUePfXTp9KrP9nU2cyFOWJENDj8cOcTs+OvvzNZpWXnnxtbdZKdH5Po8PdH/+KBrf/zAgAAAAAcNDseVEx9OsMOnzxOlEtG6fLyDKj85CdHT700miZ3+uxrY2/+OjX6w6FmRJE6OfpXQ+S3xX9xjtqBxFrSJ198qWd5vNH28HM/7xl//d1r97UihxwR5Sff/PLkz0ef899r8sPpET+QaPPDldOjIyk/Qnj7l7Qyrpg6m7kwP/jd0eeO+/vPUMJ/wp0++86kePTU6FNpf7d/fWO894WR1Sc59e7bU/rwC6MjKf8I/3JmPPniQU4PAwAAAADYlB1Of3Knr87R4LEhIkofPxpZkQGVHP6B31LXhn789KD88pPxRaJLH12sDz7bauWnRx4b8GY+3rAGekuSw6f94/P00GBUlRdX1G3nFoqMx3v99v3AyF+Pjq4cqfA/17OtGOahJx+OkPIff0oPdux/4muDqlSYX30O7jQbHn2pPajy4P2HVDmf35nPBwAAAACwn3Z4pMLPffKzevqOHYlMXlqVAdXUk4hSkfwGfexob3s770lF6NLiAtEaCUi7JD3y2MD02PlM5vytcp/mFyssel9yrRcPDQ3RsgwoLgbX2E0beuiB5YleRId36gMAAAAAAOyjnQ0qpj6dYR75rfMmdus5oFQ5n1fz+QNQzKwN/Xh0yE/Kmvv925nfryp4WKfk2n8Viw772VMXz73yu7V39Esy4n6ilDt99rWxHf4UAAAAAAD7YieDCj9HqLM5Pjv++juT1z/Njazu+18olMXgkw+w443Bf/tdYZ6oVbS9sFih2Nd6dvC8NssPLfzW/9VpeqCj7uLB+w+9P5bPt89yydSvxmaTwy+8sFFxxMVzY9no10+jiAIAAAAA7jg7WVNx+fOvxOCxzg7+vmNHIlS6drVZw5Cf/OAzIiJyp8++P8MP33+ciB568uHgzPnWyta58U9mxeDjfhDSk4h6M1c7X9L5drI4f6s17bZUriBzF/4+k2mvUidzC1WKpZZHNdrQsX6aOd/aJ3v1eoUYERGlk1FqF1HI3IU/3OIkexLRjrPKjX8yu+kTBAAAAAA40MTLL7+8IweSuQu//XCh/9EffG2pPoJY+Cif++jKp6XE4+bMxdgPnjPGXv/NBxMTn1yrHjk56hdPs/DRh4/UPvrNbz6YmJiYmM73fbc9p1Pk3kPVK//now8+npiY+L/ekz86krtSTJ545N4gUW+ifuUPU7+fmPjcuu/he4OdZxK5V5uZmPzjxMe/n4s+fix4bfJSNfFwa59ax4+5yxPXtaPfOHb06CNPBL/8138bm/h4YmLi42n7xKkXVw2tpI8fqV/+8D8+mJiYmJjI9h2PZ4vRhx45HA7fO7C0/ZMbA08dr16fM+9bdZLfHLm//VkmJm4e+tax2rVs8Mgjh8M78s8PAAAAALBvmFJqv88BAAAAAABuYztWU3Hz5s333ntvp4529zh9+nQqhToLAAAAALiNYaQCAAAAAAC6ssOL3wEAAAAAwN0GQQUAAAAAAHQFQQUAAAAAAHQFQQUAAAAAAHQFQQUAAAAAAHQFQQUAAAAAAHQFQQUAAAAAAHQFQQUAAAAAAHQFQQUAAAAAAHTl/wOJYp+Wvn/GBwAAAABJRU5ErkJggg==)
faster than they can understand this…
| state | timestamp | cases | total_population |
|---|---|---|---|
| AK | 2022-01-25 04:00:00 | 203110 | 731545 |
| AL | 2022-01-25 04:00:00 | 1153149 | 4903185 |
| AR | 2022-01-25 04:00:00 | 738652 | 3017804 |
| AZ | 2022-01-25 04:00:00 | 1767303 | 7278717 |
| CA | 2022-01-25 04:00:00 | 7862003 | 39512223 |
| CO | 2022-01-25 04:00:00 | 1207991 | 5758736 |
| CT | 2022-01-25 04:00:00 | 683731 | 3565287 |
| DC | 2022-01-25 04:00:00 | 128550 | 705749 |
| DE | 2022-01-25 04:00:00 | 241367 | 973764 |
| FL | 2022-01-25 04:00:00 | 5324438 | 21477737 |
| GA | 2022-01-25 04:00:00 | 2202292 | 10617423 |
| HI | 2022-01-25 04:00:00 | 204933 | 1415872 |
Example 2
This table shows average monthly revenue for Acme products.
| category | product | Jan | Feb | Mar | Apr | May | Jun | Jul | Aug | Sep | Oct | Nov | Dec |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| party supplies | balloons | 892 | 1557 | 1320 | 972 | 1309 | 1174 | 1153 | 1138 | 1275 | 1178 | 1325 | 1422 |
| party supplies | confetti | 1271 | 1311 | 829 | 1020 | 1233 | 1061 | 1088 | 1395 | 1376 | 1152 | 1568 | 1412 |
| party supplies | party hats | 1338 | 1497 | 1445 | 956 | 1372 | 1482 | 1048 | 877 | 1404 | 1030 | 1458 | 1547 |
| party supplies | wrapping paper | 1396 | 1026 | 932 | 891 | 1364 | 896 | 900 | 1221 | 1146 | 967 | 1394 | 1507 |
| school supplies | backpacks | 1802 | 1773 | 1611 | 1723 | 1799 | 1730 | 1813 | 1676 | 1748 | 1652 | 1819 | 1759 |
| school supplies | notebooks | 1153 | 1471 | 1541 | 1371 | 1592 | 1514 | 1725 | 1702 | 1457 | 1604 | 1729 | 1279 |
| school supplies | pencils | 1679 | 1304 | 1054 | 1259 | 1425 | 1608 | 1972 | 1811 | 1610 | 1004 | 1417 | 1283 |
| school supplies | staplers | 1074 | 1708 | 1439 | 1154 | 1551 | 1099 | 1793 | 1601 | 1647 | 1666 | 1389 | 1511 |
Use the table above to answer these questions:
- What product and month had the highest average monthly revenue?
- What product and month had the lowest average monthly revenue?
Table as Heat Map
Let’s display the table as a heat map, with larger numbers represented by darker color cells. Now how quickly can we answer those same two questions? What patterns can we see in the heat map that were not obvious in the table above?
Example 3
In 1973, Francis Anscombe published “Graphs in statistical analysis”, a paper describing four bivariate datasets with identical means, variances, and correlations.
| x1 | y1 | x2 | y2 | x3 | y3 | x4 | y4 |
|---|---|---|---|---|---|---|---|
| 10 | 8.04 | 10 | 9.14 | 10 | 7.46 | 8 | 6.58 |
| 8 | 6.95 | 8 | 8.14 | 8 | 6.77 | 8 | 5.76 |
| 13 | 7.58 | 13 | 8.74 | 13 | 12.74 | 8 | 7.71 |
| 9 | 8.81 | 9 | 8.77 | 9 | 7.11 | 8 | 8.84 |
| 11 | 8.33 | 11 | 9.26 | 11 | 7.81 | 8 | 8.47 |
| 14 | 9.96 | 14 | 8.10 | 14 | 8.84 | 8 | 7.04 |
| 6 | 7.24 | 6 | 6.13 | 6 | 6.08 | 8 | 5.25 |
| 4 | 4.26 | 4 | 3.10 | 4 | 5.39 | 19 | 12.50 |
| 12 | 10.84 | 12 | 9.13 | 12 | 8.15 | 8 | 5.56 |
| 7 | 4.82 | 7 | 7.26 | 7 | 6.42 | 8 | 7.91 |
| 5 | 5.68 | 5 | 4.74 | 5 | 5.73 | 8 | 6.89 |
| x1 | y1 | x2 | y2 | x3 | y3 | x4 | y4 | |
|---|---|---|---|---|---|---|---|---|
| mean | 9 | 7.50 | 9 | 7.50 | 9 | 7.50 | 9 | 7.50 |
| var | 11 | 4.13 | 11 | 4.13 | 11 | 4.12 | 11 | 4.12 |
Anscombe data as plots
Despite their identical statistics, when we actually plot the data we see the four datasets are very different.
)
Anscombe’s point? TO UNDERSTAND THE DATA, WE MUST PLOT THE DATA!
Rules for Effective Data Visualizations
Anscombe used clever thinking and simple plots to demonstrate the importance of data visualizations. But it’s not enough to just plot the data. To have impact, a plot must convey a clear message. How can we do this?
Have a purpose. Every data visualization should tell a story. What story does the Covid plot tell? What message does Anscombe’s quartet of plots convey?
Consider your audience. Avoid scientific names, acronyms, or jargon unless your audience is well-versed in that language. Use color-blind friendly colors.
Use an appropriate visualization. For example:
- Line graphs work well for showing changes in continuous data over time
- Bar charts compare counts or proportions in categorical data (pie charts get a bad rap in the data viz world, but can be useful in certain situations)
- For statistics (e.g., means) with confidence intervals, point plots with error bars are preferred over bar charts (nice explanation here)
- Scatterplots are useful for showing the relationship (correlation) between two continuous variables.
- Matrix heat maps can efficiently compare the magnitude of numbers when we have lots of data structured in table format, especially when colors have a clear connection to the numbers (e.g., scorecard data)
- Box plots, violin plots, and histograms show distributions and outliers for continuous data. Dot plots are a useful alternative when sample sizes are small.
KEEP IT SIMPLE!!!
- Every plot element and aesthetic should have a purpose.
- Avoid 3D charts unless you have good reason otherwise.
- Don’t try to cram everything into one plot (e.g., juxtapose two plots instead of adding a secondary y-axis. Nice explanation here).
Use informative text and arrows wisely. Clear and simple titles, subtitles, axes titles/labels, and annotations help convey the message of a plot. Use lines and arrows (sparingly but effectively) to emphasize important thresholds, data points or other plot features. Fonts should be large enough with good contrast (against the background) and sufficient white space to be easily readable.
Let’s Start Coding
For today’s data visualizations we will use the park visits data you’re already familiar with from the Data Wrangling module.
From subsets of the park visits and gas price data, we will create plots to:
- compare the distribution of park visits among years and unit types
- show temporal trends and unusual years in park visit patterns
- show temporal trends in ranking of the top 5 visited park units
- show the relationship between gas price and park visits
First steps:
- Load packages for today’s module
- Import and view the park visits and gas price data
- Extract a subset of the park visits data to work with
- Merge in the gas price data
- Check and clean/tidy the data
Step 1. Load packages
# This is a useful chunk of code that installs any missing packages and then loads all of them. Useful if we're sharing code with others who may have older versions of RStudio installed (more recent versions have pop-up prompts that ask if we want to install missing packages)
pkgs <- c("tidyverse", # includes dplyr, ggplot2, and other really useful packages
"magrittr", # for the `%<>%` function to reassign piping outcomes to the original variable
"viridis", # popular colorblind-friendly palette
"scales", # show_col() and scale_x_date()
"knitr") # format tables for display
installed_pkgs <- pkgs %in% installed.packages() # check which packages are not yet installed
if (length(pkgs[!installed_pkgs]) > 0) install.packages(pkgs[!installed_pkgs],dep=TRUE) # if some packages are not yet installed, go ahead and install them...
lapply(pkgs, library, character.only = TRUE) # ...then load all the packages and their dependenciesStep 2. Import and view the data
park_visits <- read_csv("https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2019/2019-09-17/national_parks.csv") # visits to US National Parks
gas_price <- read_csv("https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2019/2019-09-17/gas_price.csv") # gas prices over space and timeView() the data to get an initial sense of what we’re working with.
# Examine the data to understand data structure, data types, and potential problems
View(park_visits); summary(park_visits)
View(gas_price); summary(gas_price)Step 3. Extract a subset of data
Before we run more detailed data checks, let’s extract a subset of the park_visits data to work with.
CHALLENGE: Extract a subset of the park_visits data
Create a new data frame that contains a subset of the park_visits data. Call the new data frame park_sub. We only want to include park_visits data that meet these criteria:
- region is IM, PW or SE
- unit_type is National Monument, National Historic Site, or National Park
- year is in the range 2000 - 2015, inclusive
We also only need to keep a subset of the data columns from park_visits. The new data frame should have these columns (in this order): region, unit_type, unit_name, unit_code, year, visitors
Quick Tip!
Even though park_visits$year is classified as a character data type, R can sort these numbers-as-characters in numeric order. Try running sort(unique(park_visits$year)) in the console to see for yourself. This means we can filter park_visits for a subset of years "2000":"2015" without first converting park_visits$year to a numeric data type (actually we can even omit the quotes and simply filter for years in 2000:2015).
The new data frame should look like this when we type glimpse(park_sub) in the console:
## Rows: 1,919
## Columns: 6
## $ region <chr> "PW", "SE", "IM", "IM", "PW", "IM", "PW", "PW", "PW", "IM", ~
## $ unit_type <chr> "National Historic Site", "National Historic Site", "Nationa~
## $ unit_name <chr> "Manzanar National Historic Site", "Tuskegee Airmen National~
## $ unit_code <chr> "MANZ", "TUAI", "SAND", "WABA", "CECH", "WACO", "NPSA", "CRL~
## $ year <chr> "2000", "2003", "2010", "2002", "2013", "2015", "2002", "201~
## $ visitors <dbl> 38010, 11085, 4063, 15374, 8157, 20551, 1938, 614712, 326376~
Step 4. Merge the gas_constant data with the park_sub data frame
CHALLENGE: Join the gas_constant data (from gas_price data frame) to park_sub, matching records by year. Call the new data frame joined_dat.
The new data frame should show the additional data column gas_constant when we type glimpse(joined_dat) in the console:
## Rows: 1,919
## Columns: 7
## $ region <chr> "PW", "SE", "IM", "IM", "PW", "IM", "PW", "PW", "PW", "IM~
## $ unit_type <chr> "National Historic Site", "National Historic Site", "Nati~
## $ unit_name <chr> "Manzanar National Historic Site", "Tuskegee Airmen Natio~
## $ unit_code <chr> "MANZ", "TUAI", "SAND", "WABA", "CECH", "WACO", "NPSA", "~
## $ year <dbl> 2000, 2003, 2010, 2002, 2013, 2015, 2002, 2015, 2015, 201~
## $ visitors <dbl> 38010, 11085, 4063, 15374, 8157, 20551, 1938, 614712, 326~
## $ gas_constant <dbl> 2.02, 2.01, 3.02, 1.75, 3.62, 2.45, 1.75, 2.45, 2.45, 2.4~We now have all the data we need in a single data frame.
Step 5. Check and clean/tidy the data
We are primarily using dplyr to wrangle data and ggplot2 to plot data. These two packages are part of the Tidyverse of packages for data science. Tidyverse packages share a common design philosophy centered on ideas like breaking down complex tasks into a series of simpler functions executed in sequence (e.g., piping %>%). Functions from tidyverse packages work together most seamlessly and efficiently when data are in a tidy format. Remember that a ‘tidy’ format means data have one variable per column and one observation per row.
POLL: Which of these is ‘tidy’ data?
- Plants that Grow
- Political Twitter Likes
- Test Grades and Lunch
- A & C
- None of the above
| plant_id | age_days | height_cm | num_leaves |
|---|---|---|---|
| 1 | 160 | 21 | 11 |
| 2 | 160 | 19 | 12 |
| 3 | 160 | 32 | 16 |
| 4 | 160 | 23 | 8 |
| 1 | 190 | 24 | 12 |
| 2 | 190 | 20 | 14 |
| political_party | Cat Food Breath | Doug the Pug | Jon Pigeon |
|---|---|---|---|
| Democrat | 5201 | 4760 | 1200 |
| Republican | 6912 | 3448 | 800 |
| Independent | 3722 | 2666 | 5128 |
| student | before_lunch | after_lunch |
|---|---|---|
| Robin | 89 | 71 |
| Marl | 77 | 76 |
| Randi | 82 | 43 |
What about our joined_dat data frame? Is it in ‘tidy’ format? Yes, it is! No need to pivot_long anything today. But let’s make two final changes to the data before we begin plotting. First, we will convert region to an ordered factor with more descriptive labels. Second, we will remove park units that don’t have visitor data for every year from 2000 to 2015. We will end with a final check of the data to make sure it’s ready for plotting.
Convert region to a factor data type
Convert the variable region to a factor data type and rename and re-order the levels as Pacific West < Intermountain < and Southeast (the default order is alphabetical). This way, bar plots and color legends will order the region levels from west to east instead of alphabetically.
# NOTES:
# 1. The factor data type is similar to (and often interchangeable with) the character data type, but limits entries to pre-specified values and allows us to specify an order to those values (other than the default alphabetical order)
# 2. Factors are especially useful for ordered categories (e.g., low < medium < high) with a small number of possible values
# 3. Under the hood, `R` store factors as integer vectors representing the order of the value (with character string labels)
# 4. I usually classify categorical variables as character data types, unless I want to limit entries to a few pre-specified values, or if I want to enforce a particular order to the category levels
unique(joined_dat$region) # check what values are in `region`## [1] "PW" "SE" "IM"joined_dat$region <-
factor(joined_dat$region,
levels = c("PW", "IM", "SE"), # specifies the factor levels and an order for the levels. If we omit this argument, the factor levels will be the unique set of values in the column, ordered alphabetically
labels = c("Pacific West", "Intermountain", "Southeast")) # relabels the levels with more descriptive names
# Check the data
glimpse(joined_dat) # `region` is now a factor data type## Rows: 1,919
## Columns: 7
## $ region <fct> Pacific West, Southeast, Intermountain, Intermountain, Pa~
## $ unit_type <chr> "National Historic Site", "National Historic Site", "Nati~
## $ unit_name <chr> "Manzanar National Historic Site", "Tuskegee Airmen Natio~
## $ unit_code <chr> "MANZ", "TUAI", "SAND", "WABA", "CECH", "WACO", "NPSA", "~
## $ year <dbl> 2000, 2003, 2010, 2002, 2013, 2015, 2002, 2015, 2015, 201~
## $ visitors <dbl> 38010, 11085, 4063, 15374, 8157, 20551, 1938, 614712, 326~
## $ gas_constant <dbl> 2.02, 2.01, 3.02, 1.75, 3.62, 2.45, 1.75, 2.45, 2.45, 2.4~levels(joined_dat$region) # shows the (renamed) factor levels and their order## [1] "Pacific West" "Intermountain" "Southeast"str(joined_dat$region) # shows the values are coded as integers indicating factor level order## Factor w/ 3 levels "Pacific West",..: 1 3 2 2 1 2 1 1 1 2 ...summary(joined_dat) # the summary for `region` now counts the number of times each factor level occurs## region unit_type unit_name unit_code
## Pacific West :527 Length:1919 Length:1919 Length:1919
## Intermountain:965 Class :character Class :character Class :character
## Southeast :427 Mode :character Mode :character Mode :character
##
##
##
## year visitors gas_constant
## Min. :2000 Min. : 0 Min. :1.750
## 1st Qu.:2004 1st Qu.: 57344 1st Qu.:2.320
## Median :2008 Median : 187397 Median :3.000
## Mean :2008 Mean : 602483 Mean :2.828
## 3rd Qu.:2012 3rd Qu.: 652166 3rd Qu.:3.610
## Max. :2015 Max. :10712674 Max. :3.800Remove park units with missing years of data
We will be creating plots to see if trends in park visits differ by region or unit type, so we want to exclude park units for which we don’t have complete data for the years 2000 to 2015.
Finding which park units have missing years of data is not as simple as filtering for NA’s in joined_dat because missing park unit-years are not represented by records with NA for the visitor count–the missing data are simply not in the data frame.
We could use tidyr’s expand_grid function (and then join the template data frame with joined_dat) to create the missing records with an NA entry for visitors. We would then delete any park units that have at least one NA in the visitors column. expand_grid is a useful function to know because we sometimes need to add in missing records (and, for example, fill in zero values to indicate we went to the site but found zero saplings) so we can get correct summary statistics such as means and standard deviations (that should not be treating zero counts as missing data).
In our case, we have a simpler option. We can count the unique years of data for each park unit and delete park units that have less than 16 unique years of data.
(Note that we should also check for park unit-year combinations with multiple records. If these are duplicate records, the duplicates can be deleted. If these records have different visitor counts for the same park unit-year, we would have to find out which record is the “correct” one).
length(2000:2015) # if we can't figure this out in our head, we probably need more sleep. But here is proof that there are 16 years in the range 2000 to 2015, inclusive## [1] 16# Check if we have multiple records for any combination of unit_code and year (e.g., ARCH should not have multiple records for the year 2016)
# Short example to show how `duplicated` works
test <- data.frame(x = c(1, 1, 3, 2, 2), y = c(7, 2, 4, 1, 1), z = 1:5)
duplicated(test[c("x", "y")])## [1] FALSE FALSE FALSE FALSE TRUEtest[duplicated(test[c("x", "y")]), ]## x y z
## 5 2 1 5# Now use `duplicated` on our data
joined_dat[duplicated(joined_dat[c("unit_code", "year")]), ] # good, no duplicates## # A tibble: 0 x 7
## # ... with 7 variables: region <fct>, unit_type <chr>, unit_name <chr>,
## # unit_code <chr>, year <dbl>, visitors <dbl>, gas_constant <dbl># Identify park units with less than 16 years of data
remove_units <- joined_dat %>%
count(unit_code) %>% # count() combines group_by() and summarize(n = n()) in a single function.
dplyr::filter(n != 16) # sort the data frame from lowest to highest count of unique years
# View(joined_dat) to check results before you remove these park units from the data
# Remove those park units from the data
viz_dat <- joined_dat %>%
dplyr::filter(!unit_code %in% remove_units$unit_code)
# View(viz_dat) to make sure the park units have been removed. The final data frame should have 1856 rows and 7 columns.# If you want to save `viz_dat` to your current working directory, run the code below. This code saves `viz_dat` as an .RDS file, which just means it saves it as an `R` object instead of as a .csv file.
# When we save a data frame as an `R` object we retain the accompanying attribute information (e.g., the order of factor levels for `region`)--this information would not be saved in a .csv file.
saveRDS(viz_dat, "viz_dat.RDS")
# To read the RDS file from your working directory and assign it to `viz_dat`, use the code below. You can assign it to any name you want. If you're reading it from a location other than your working current directory, provide the file path in the function argument.
viz_dat <- readRDS("viz_dat.RDS")Note that a report of our findings should explain how and why we selected this subset of data to work with.
Final data check
Always do a final check of the data to look for missing data, duplicated records, funny numbers or category levels, etc. It helps to create summary tables that decompose the data in various ways. These tables can identify oddities (e.g., category levels with only one record) and also give us an idea of the sample sizes we’ll be working with if we create plots with subsets of the data.
CHALLENGE: Create a table (data frame) that shows the number of park units in each combination of region and. Extra credit if you can reformat the table to wide instead of long format.
The new data frame (if you have it in wide format) should look like this:
## # A tibble: 3 x 4
## unit_type `Pacific West` Intermountain Southeast
## <chr> <int> <int> <int>
## 1 National Historic Site 7 9 12
## 2 National Monument 10 35 8
## 3 National Park 17 18 7
We see from this table that we have data for 34 National Monuments in the Intermountain region and that the smallest number of park units in any region-type is 7. We’re ready to start plotting.
Plotting with Base R
# This tab requires the `viz_dat` data frame created earlier in the module. If you don't have this data frames in the global environment, run the code below.
viz_dat <- readRDS("viz_dat.RDS") # this file is provided with the training materials. If this file is not in your current working directory, provide the file path in the function argumentPlots can reveal data issues that we had overlooked with our data checks. I sometimes use base R plotting functions to create quick-and-dirty plots for this purpose. If I want anything more than a very basic plot, I find ggplot functions easier to remember and use.
In addition to quick plots, base R plotting functions have another important use. Many R packages have functions that create specific data visualizations when R’s generic plot() is applied to an object of a particular class. For example, lm() is a statistical function that fits a general linear model to data and saves the model as an object of class lm. Run the code below to see what happens when we apply plot() to an object of class lm.
# NOTE: `lm` is part of the `stats` package that is automatically installed with `R`. This code comes from the Examples section of `?lm`.
# Create fake data
ctl <- c(4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14)
trt <- c(4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69)
group <- gl(2, 10, 20, labels = c("Ctl","Trt"))
weight <- c(ctl, trt)
# Simple linear regression of `weight` as a function of group
mod <- lm(weight ~ group) # assign model results to the variable `mod`
summary(mod) # model summary
class(mod) # class of `mod` is `lm`
plot(mod) # residual plotsWe aren’t going to discuss statistics or residual plots today, but it’s useful to know these common uses for plot().
Feeding Data to Base R Plotting Functions
Look in the help file for the R graphics package that is automatically installed with R (in RStudio’s ‘Packages’ tab, click the hyperlink for the package called ‘graphics’). Here you will find many functions for base R graphics. For example, the boxplot() function draws box plots.
Typically, we can feed data to most base R plotting functions in one of two ways:
Provide the data as a formula of the form
y ~ xwhereyis a vector of data for the y-axis andxis a vector of data for the x-axis. For example,plot(y ~ x, data = my_data)orplot(my_data$y ~ my_data$x).Provide the vectors as separate arguments. For example,
plot(x = my_data$x, y = my_data$y)When providing vectors as separate arguments, each vector must be a standalone vector. That is, we can’t just give column names and provide the data frame via a data argument.
With the non-formula approach, the argument names may not be x and y. For example, with bar plots, the argument name height is used instead of y. Refer to the help (e.g., ?barplot) for each plotting function to see what arguments it accepts.
I tend to use the formula approach for feeding data to base R plotting functions, because it is more consistent across the different types of plots. In other words, I’m too lazy to look up what argument names to use with each plot type for the non-formula approach.
Example: Box Plots
Let’s use base R functions to create a box plot of park visitors by year, for a subset of viz_dat that meets these criteria:
- region is Intermountain
- unit_type is National Monument
- year is 2000, 2005, 2010, or 2015 (each year will be a separate box plot)
We want to know how much spread there is (across park units) in the number of visitors to National Monuments of the Intermountain region. We also want to know if visitor counts differed much among the four years.
# Check out `?boxplot`
# Subset the data
subdat <- viz_dat %>%
dplyr::filter(
year %in% c(2000, 2005, 2010, 2015),
unit_type == "National Monument",
region == "Intermountain")
boxplot(visitors ~ year, data = subdat) # the basic boxplot, with data provided as a formula
# An alternative way to feed the data as a formula
# boxplot(subdat$visitors ~ subdat$year) # no need for `data` argumentThis plot is okay, but even for a basic plot I would want to make the y-axis labels more understandable. Because the visitor numbers are large, R is using scientific notation for the y-axis labels. So instead of 200,000 we see 2e+05. Let’s change that in the box plot and add plot and axes titles.
# Save this boxplot to `p_box` so we can look at the underlying structure
p_box <- boxplot(visitors/1000 ~ year, # so now, 200,000 will be shown as 200 on the y-axis
data = subdat,
main = "National Monuments in the Intermountain Region", # plot title
xlab = "Year", # x-axis title
ylab = "Number of visitors, in 1000's") # y-axis title
Not bad for basic box plots. Now look at the data underlying the box plots. In ?boxplot, the section titled ‘Value’ explains what these numbers mean. For any
p_box## $stats
## [,1] [,2] [,3] [,4]
## [1,] 3.131 2.919 4.3500 9.492
## [2,] 61.170 53.521 52.5660 53.165
## [3,] 112.611 103.570 103.8875 94.931
## [4,] 258.306 280.068 221.0830 222.723
## [5,] 550.657 505.158 470.9210 416.635
##
## $n
## [1] 34 34 34 34
##
## $conf
## [,1] [,2] [,3] [,4]
## [1,] 59.1935 42.18307 58.22483 48.98625
## [2,] 166.0285 164.95693 149.55017 140.87575
##
## $out
## [1] 789.131 848.348 622.320 830.253 525.831 578.554 827.247 793.601 497.506
## [10] 478.833 813.686 588.006
##
## $group
## [1] 1 1 2 2 3 3 3 4 4 4 4 4
##
## $names
## [1] "2000" "2005" "2010" "2015"From the boxplots and numbers, it looks like visitation was generally down in 2010 and 2015 compared to 2000 and 2005. The spread of visitor counts across these 34 National Monuments is huge, though.
The code below extracts the minimum and maximum visitor counts in subdat for the year 2000.
range(subdat$visitors[subdat$year==2000]) ## [1] 3131 848348For the Intermountain region in the year 2000, the number of visitors (N = 848348) to the most-visited National Monument was 270 times the number of visitors (N = 3131) to the least-visited National Monument!
Modifying Base R Plots
Type ?par in the console to see the many graphical parameters (arguments) that can be used to modify base R plots. Unless you plan to use base R functions to create complex plots, you can probably get by with remembering a small handful of useful graphical parameters, such as pch = to change scatterplot point types, col = to set color, and a variety of cex arguments to change font size and point size.
Let’s create a basic scatterplot and use some of these graphical parameters to gussy up the plot.
Example: Scatterplots
CHALLENGE: Create a scatterplot showing the relationship between park visitor count (y-axis) and gas price (x-axis) for Arches National Park (unit code is ARCH).
We want to know if annual visitation to Arches National Park is correlated with gas price (annual national average).
The data for this plot should be the 16 records of data for Arches National Park (one data point per year). Assign this subset of data to the name ‘arch_dat’ because we will be using it again for later plots.
Make these modifications to the default scatterplot:
Show visitor count as 1000’s of visitors (i.e., divide the visitor count by 1000) as we did for the box plots.
Use the graphical parameter
pchto change the point type from the default open circle to a solid circle (see?pointsto find the number corresponding to a filled circle. The examples at the end of?pointsshow how to usepchas a plot argument)Increase point size to 1.5 instead of the default of 1
The final plot should show this relationship: 
This plot shows that annual visitor counts increased as gas prices increased. That’s not the relationship we might have expected to see. But correlation is not causation, and the resolution of data influences the pattern we see. Maybe if we had used average June - July gas prices within 330 km of Arches National Park (instead of annual national averages) we might have seen a different relationship with visitor counts. Maybe, maybe not.
The scatterplot shows another interesting bit of information. What is that funny point that is marching to its own beat?! Gas price was low (~2.50 per gallon in 2015 dollars) and the visitor count was really high, ~ 1.4 million. A quick look at the data shows that this funny point was for the year 2015.
It would be interesting to see if 2015 was a high visitation year for other park units as well. We will save that question for ggplot.
For now, let’s see if we gain any additional insight from plots showing trends in visitor counts at Arches National Park and trends in national gas prices. In the plots below, the data point for year 2015 is shown in red. To do that, I added this argument to the plot function: col = ifelse(arch_dat$year == 2015, "red", "black")
The plot above shows that from 2004 to 2015, visitor numbers at Arches National Park increased every year. Visitation REALLY jumped in 2014 and 2015.
The plot above shows that gas prices have gone up and down between 2000 and 2015. There was a huge drop in gas price between 2008 and 2009 and again between 2014 and 2015.
These scatterplots give interesting insight on trends in national gas prices and trends in visitation at Arches National Park. The pattern we saw in our initial scatterplot–positive correlation between annual visitation and national average gas price–was almost certainly a spurious correlation. It’s important to look at data from many different directions and to not jump to conclusions based on a narrow exploration of the data.Example: Line Graph
Line graphs are often used for showing changes in continuous data over time. Let’s make one last plot in base R before we move on to ggplot.
CHALLENGE: Re-create the plot of trends in visitation at Arches National Park, using a line graph instead of a scatterplot. If you’re feeling ambitious, make it a line graph with overplotted points.
The final plot (with overplotted points) should show this relationship: 
Building a Basic ggplot()
# This tab requires the `viz_dat` and `arch_dat` data frames created earlier in the module. If you don't have these data frames in the global environment, run the code below.
viz_dat <- readRDS("viz_dat.RDS") # this file is provided with the training materials. If this file is not in your current working directory, provide the file path in the function argument
arch_dat <- subset(viz_dat, unit_code == "ARCH") # subset of data that is for Arches National parkIf you search the Internet for information on plotting in R, you will quickly learn that ggplot is the most popular R package for plotting. It takes a little effort to learn how the pieces of a ggplot() object fit together. But once you get the hang of it, you will be able to create a large variety of attractive plots with just a few lines of R code.
Reference documents and cheatsheets can be very helpful while you are learning ggplot().
We will build our first ggplot() step-by-step to demonstrate how each component contributes to the final plot. Our first ggplot() will re-create the line graph of trends in visitation at Arches National Park.
Step 1. Create a “Master” Template for the Plot
- identify the data (must be a data frame or tibble) [REQUIRED]
- set default aesthetic mappings
An aesthetic mapping describes how variables in the data frame should be represented in the plot. Any column of data that will be used to create the plot must be mapped to an aesthetic. Examples of aesthetics include x, y, color, fill, point size, point shape, and line type.
From the ‘arch_dat’ data we will map ‘visitors’ (in 1000’s) to the y variable and ‘year’ to the x variable. This will set up default axes and axes titles for the plot.
Master template with default aesthetics
# NOTES:
# 1. The first argument provided to `ggplot()` is assumed to be the data frame, so it's not necessary to include the argument name `data =`.
# 2. We are assigning the plot to a variable `p` so we can build on this `ggplot()` master template in the next step.
# 3. The parentheses around the line of code is a shortcut for `print`. Without it, the plot would assign to `p` but not print in the plots pane.
(p <- ggplot(data = arch_dat, aes(x = year, y = visitors/1000)))
summary(p) # summary of the information contained in the plot object## data: region, unit_type, unit_name, unit_code, year, visitors,
## gas_constant [16x7]
## mapping: x = ~year, y = ~visitors/1000
## faceting: <ggproto object: Class FacetNull, Facet, gg>
## compute_layout: function
## draw_back: function
## draw_front: function
## draw_labels: function
## draw_panels: function
## finish_data: function
## init_scales: function
## map_data: function
## params: list
## setup_data: function
## setup_params: function
## shrink: TRUE
## train_scales: function
## vars: function
## super: <ggproto object: Class FacetNull, Facet, gg>p$data # the plot object is self-contained. The underlying data frame is saved a list element in the plot object. ## # A tibble: 16 x 7
## region unit_type unit_name unit_code year visitors gas_constant
## <fct> <chr> <chr> <chr> <dbl> <dbl> <dbl>
## 1 Intermountain National Park Arches Nat~ ARCH 2015 1399247 2.45
## 2 Intermountain National Park Arches Nat~ ARCH 2014 1284767 3.4
## 3 Intermountain National Park Arches Nat~ ARCH 2013 1082866 3.62
## 4 Intermountain National Park Arches Nat~ ARCH 2012 1070577 3.8
## 5 Intermountain National Park Arches Nat~ ARCH 2011 1040758 3.75
## 6 Intermountain National Park Arches Nat~ ARCH 2010 1014405 3.02
## 7 Intermountain National Park Arches Nat~ ARCH 2009 996312 2.58
## 8 Intermountain National Park Arches Nat~ ARCH 2008 928795 3.61
## 9 Intermountain National Park Arches Nat~ ARCH 2007 860181 3.16
## 10 Intermountain National Park Arches Nat~ ARCH 2006 833049 3
## 11 Intermountain National Park Arches Nat~ ARCH 2005 781670 2.74
## 12 Intermountain National Park Arches Nat~ ARCH 2004 733131 2.32
## 13 Intermountain National Park Arches Nat~ ARCH 2003 757781 2.01
## 14 Intermountain National Park Arches Nat~ ARCH 2002 769672 1.75
## 15 Intermountain National Park Arches Nat~ ARCH 2001 754026 1.91
## 16 Intermountain National Park Arches Nat~ ARCH 2000 786429 2.02Step 2. Add Geometry Layer(s)
The geometry layer determines how the data will be represented in the plot. We can generally think of this as the plot type. For example, geom_hist() creates a histogram and geom_line() creates a line graph.
If we had set default aesthetic mappings, the geometry layer will adopt those mappings unless additional (overriding) aesthetics are set in the layer itself. If we had not set default aesthetic mappings, they will need to be set within the layer. At a minimum we need to specify which data frame columns hold the x and (for most plot types) y aesthetics. Except for these two aesthetics, ggplot() will use its default values for any aesthetics not defined by the user.
Each geometry layer has its own set of accepted aesthetics. For example, the aesthetic linetype cannot be used in a geom_point() (point plot) layer. We cover plot aesthetics in more detail in the tab titled ‘Customizing a ggplot()’.
Add a line graph geometry layer
# NOTE: We combine `ggplot()` layers with `+`, not with `%>%`. With a `ggplot()` object we are not piping sequential results but rather layering pieces. In this example we are building our `ggplot()` object in a particular order, but we could actually reorder the pieces in any way (as long as the master template is set first). The only difference would be that if different layers have conflicting information, information in later layers overrides information in earlier layers.
p + geom_line()
Example of aesthetics set in the geometry layer
If we had NOT set the x and y aesthetics in the plot’s master template we could have set them within the geometry layer like this:
# Example of setting aesthetic mappings within the layer only
ggplot(data = arch_dat) +
geom_line(aes(x = year, y = visitors/1000))Add multiple geometry layers
We can add include multiple geometry layers in a single ggplot() object. For example, to overlay a point plot on the line graph we would add a geom_point() layer. If we had set default aesthetic mappings, this additional geometry would adopt those mappings unless additional (overriding) aesthetics were set in the geom_point() layer itself.
p + geom_line() + geom_point()
If we had only set aesthetic mappings within the geom_line() layer, we would also need to provide geom_point() with its own aesthetic mappings – there are no default aesthetics for it to adopt.
The code below generates an error message: Error: geom_point requires the following missing aesthetics: x and y
# Example of setting aesthetic mappings within the layer only
ggplot(data = arch_dat) +
geom_line(aes(x = year, y = visitors/1000)) +
geom_point()
Geometry layers are shortcuts for the layer() function
Geometry layers can alternatively be specified with the generic plot layer() function. Geometry layers are basically the layer() function with some pre-defined (default) arguments. For example, geom_point() is the same as layer(geom = "point", stat = "identity", position = "identity"). See ?layer for more information.
p + layer(geom = "line", stat = "identity", position = "identity") # this creates the same line graph we had created using `geom_line()`
Step 3. Use scale_...() Functions to Fine-Tune Aesthetics
All aesthetic mappings can be fine-tuned with a scale_...() function. For example, the aesthetic mapping for x can be modified with scale_x_continuous() if the data are continuous, scale_x_discrete() if the data are categorical, and scale_x_date() if the data class is data/time. ggplot provides default scale arguments for all aesthetic mappings. If we don’t explicitly define a scale for an aesthetic, ggplot will use its default scale settings.
Modify the y-axis with scale_y_continous()
# NOTES:
# 1. Our plot code is getting long, so we will separate out the components to improve readability
# 2. We will assign this plot object to a new variable, `p2` so we can just add to `p2` in the next step instead of re-typing all these lines
# 3. `ggplot()` also provides shortcuts for some scaling arguments. For example, if we just want to set y-axis limits we can add the layer `ylim (600, 1500)` instead of setting limits via `scale_y_continuous()`
# Check out `?scale_y_continuous` for ways to modify the y-axis
(p2 <- p +
geom_line() +
geom_point() +
scale_y_continuous(
name = "Number of visitors, in 1000's", # y-axis title
limits = c(600, 1500), # minimum and maximum values for y-axis
breaks = seq(600, 1500, by = 100)) # label the axis at 600, 700, 800... up to 1500
)
Step 4. Use labs() to Specify Plot Labels
This useful layer allows us to set various plot labels, such as the plot title and subtitle, x- and y-axis titles, and alternative text and plot tags. Some of these labels can be set in other ways rather than via labs(). For example, we had already set the y-axis title in the first argument of scale_y_continous().
If we want to remove a default label from the plot, we can set the argument to NULL. For example, labs(x = NULL) removes the default x-axis title.
Add plot labels
(p3 <- p2 +
labs(
title = "Arches National Park Visitation, 2000 - 2015", # plot title
subtitle = "Visitor counts have increased each year since 2004", # plot subtitle
x = "Year") # x-axis title. We had already set the y-axis title with `scale_y_continous()`.
)
Step 5. Set a Complete Theme and Plot Element Themes
Use theme() to customize all non-data components of the plot. Conveniently, ggplot() provides complete themes that modify many element themes at once. The default complete theme is theme_grey. See ?theme for all the ways we can modify the non-data plot elements. This help page also provides a hyperlink to a help page for complete themes (I’m sure there is a way to get there more directly, but I couldn’t figure it out).
Explore ggplot() complete themes
p3 +
theme_minimal() # try out different themes! theme_classic(), theme_void(), theme_linedraw()... 
Create a FiveThirtyEight-inspired plot
For more fun, install and load the ggthemes package, which allows us to apply plot styles inspired by the Economist and the popular FiveThirtyEight website, among other choices.
p3 +
ggthemes::theme_fivethirtyeight()
Modify theme elements
In addition to setting complete themes, we can individually modify element themes. If the element themes layer is set AFTER the complete themes layer, the element themes will override any conflicting default arguments in the complete theme. Some of these element these can alternatively be set in a scale_...() function. For example, to move the y-axis to the right side of the plot, use scale_y_continuous(position = "right") or theme()
It’s a bit difficult to understand the help page for ?themes. I usually resort to an Internet search to figure out how to set an element theme.
p3 +
ggthemes::theme_fivethirtyeight() +
theme(
plot.subtitle = element_text(color = "red", size = 14, hjust=0.5)) # change plot subtitle font color and size, and center the subtitle
Step 6. Split a Plot into Many Subplots
We can use facet_grid() or facet_wrap() to split a ggplot() object into multiple subplots, each representing a different group level. This is a nice way to simplify data-heavy plots with a single line of code.
With facet_grid we can specify different grouping variables for a row of plots versus a column of plots. For example, we can say that each row of subplots should represent a different region and each column of subplots should represent a different unit type. With facet_wrapped, subplots are “wrapped” in order from left-to-right and top-to-bottom, just like text is wrapped on a page.
Create a faceted plot
# NOTES:
# 1. For this plot we need to use the `viz_dat` dataset so we can facet by park unit name
# 2. In the `facet_wrap` arguments, `ncol = 1` sets a single column so plots are stacked on top of each other and `scales = "free_y"` allows the y-axis range to differ across the subplots (this is useful if the range of y-axis values is very different across subplots and we are more interested in comparing trends than in comparing absolute numbers among the subplots)
subdat <- subset(viz_dat, unit_code %in% c("ARCH", "GRCA", "YOSE")) # plot data for Arches National Park (ARCH), Grand Canyon National Park (GRCA), and Yosemite National Park (YOSE)
ggplot(subdat, aes(x = year, y = visitors/1000)) +
geom_line() +
geom_point() +
labs(
title = "National Park Visitation, 2000 - 2015",
x = "Year",
y = "Number of visitors, in 1000's") +
theme_bw() +
theme(plot.title = element_text(hjust=0.5)) + # center the plot title
facet_wrap(vars(unit_name), ncol = 1, scales = "free_y")
We were wondering if 2015 was a very high visitation year for parks other than Arches National Park. It looks like it was for Grand Canyon and Yosemite National Parks too!
Summary of Useful ggplot() Components
These are not the ONLY ggplot() components, but they are the ones most useful to know (IMHO)
Data - a data frame (or tibble). Specified as the first ggplot() argument. Any aesthetics set at this level are part of the “master” template for the plot.
Aesthetic mapping - the variables represented in the plot (and how), e.g., x, y, color, fill, size, shape, line type. These aesthetics, except x & y, can also be set to a constant value by specifying outside the aes(). Aesthetics can be defined in the master template or in a geometry layer.
Geometry layer - the plot/mark layers, e.g., geom_point, geom_line, geom_boxplot. Can specify multiple such layers in a plot.
Scales - used to fine-tune aesthetics. Takes the form ‘scale_xxx_yyy’, where ‘xxx’ refers to the aesthetic and ‘yyy’ how the aesthetic values will be specified, e.g., discrete, continuous, manual. The associated legend for the aesthetic can also be modified here.
Labs - sets the plot labels, including plot title, subtitle, caption, and axis labels. Importantly, this is also where you would set alternative text and plot tags, for Section 508 Compliance.
Themes - complete themes (e.g., theme_bw, theme_classic) can be used to control all non-data display. Themes can also be used to more specifically modify aspects of the panel background, legend elements, label text size and color, etc.
Facets - subset the data by grouping variables and generate a separate subplot for each grouping variable, putting the graphs all on the same page.
Vizualizing Spatial Data
GIS in R
If you’ve tried to do GIS and make maps in R even a few years ago, you probably encountered the same frustrations I did. There were a ton of packages out there, each with its own file format and coding convention, and each package rarely had all the features I needed. It was not easy to navigate, and I often found myself just exporting my work out of R and doing the rest in ArcGIS… Enter the sf and tmap packages, which are the latest and greatest R packages devoted to GIS and map making! Most of the frustrations I had with earlier packages have been resolved with these two packages.
The sf package is the workhorse for anything you need to do with spatial vector data. File types with sf are called simple features, which follow a set of GIS standards that are becoming the universal data model for vector-based spatial data in R and that most GIS-based packages in R now employ. Simple features are also now the standard for open source GIS in general. That means if you’re trying to troubleshoot something with simple features, you’ll need to specify that it’s for R, rather than PostGIS or some other implementation. The sf package is also superseding the rgdal package, which used to be the more common data model in R and open source GIS. The more I use this package, the more impressed I am with how intuitive it is to use, and how well documented the functions are. For vector data, I have yet to need to perform a task that sf couldn’t do.
The main application of tmap package is making maps, and it was designed using a grammar of graphics philosophy similar to ggplot2. In practice for tmap, it means that maps are built as a collection of many small functions/layers that get added together, and order matters. There are also tmap-enabled functions that you can use in ggplot2 plots too, but you can do a lot more in tmap. I also prefer the look of tmap’s built-in compass, legends, etc. over the ones available in ggspatial, which is an add-on package for plotting spatial data in ggplot2.
The raster package is also an excellent package for analyzing/processing raster data. For large jobs, I find the raster package easier to work with than ESRI tools, and it tends to run a lot faster than ESRI built-in tools (I haven’t compared with python).
Finally, the leaflet package in R allows you to create interactive maps as html widgets. These are often included in R shiny apps, but they can also be used in R Markdown with HTML output (for more on that, attend next Wednesday’s session on R Markdown). An internet connection is required for the maps to be interactive. Without an internet connection the map will show the default extent defined by the code.
leaflet package is relatively easy to use for basic mapping. For more advanced features or to customize it further, you often end up having to code in JavaScript, which is the language leaflet was originally programmed in. There’s also a lot more online help available for the JavaScript version of the leaflet library than the R version. If you’re really stumped about something, you may find your answer with the JavaScript help.
Using sf and tmap
My two favorite features of sf are 1) the attribute table of a simple feature (sf’s equivalent of a shapefile) is a dataframe in R, and 2) package functions are pipeable with tidyverse functions. That means, if you want to delete points in your layer, you can use dplyr’s filter() function to filter out those points. The sf package will update the geometry of the layer to remove the points that were filtered out.
To demonstrate the use of sf and tmap, I’m going to generate a simple GRTS sample using spsurvey, which now connects to sf instead of sp and rgdal. Then we’ll filter out points that don’t fall in a forest habitat to demonstrate how to work with sf data in R. Finally we’ll plot the results using tmap.
I wasn’t able to figure out how to post a shapefile that you could easily download in R with a url. I emailed them to you as data.zip the previous week. I also posted all of the files to our Scientists Team, which can be downloaded using the following link: https://doimspp.sharepoint.com/:f:/s/ScientistsTraining2022/EiXOTYV1l4RCk1sMr5yXhZUB1ZFEaAlAN4elvsYbBfYHRg?e=ktVy5n. To follow along, you’ll need to download (and unzip if you’re using the email attachment) the shapefiles and save them to your working directory.
Step 1. Load and prep shapefiles using sf
The code chunk below has a lot going on, but bare with me. This is the easiest way I’ve found to customize the layers in tmap plots. The basic steps in the code chunk below are:- Load libraries
- Set a seed for GRTS sample. As long as you save the seed and layers used to generate the GRTS sample, you should be able to recreate the same GRTS sample.
- Read in park boundary and vegetation map layers.
- View quick plot of layer in base R, just to make sure it looks right.
- Clip the veg layer to the park boundary using intersection, so map looks better.
- Look at the attribute tables for each layer.
- Add 2 columns to the veg layer to simplify the habitats (simp_veg) and set their colors (fill_col)
Once those steps are completed, we’re ready to generate a GRTS sample and start making a map. Note that I’m using 6-digit hex colors (i.e., “#ffffff” is white) to define the fill color for each habitat type. To find your own or see what color these look like, check out htmlcolorcodes.com
#install.packages(c("sf", "spsurvey"))
library(dplyr) # for filter, select, mutate and %>%
library(sf)
library(tmap)
library(spsurvey)
# Generate a random number for the seed
sample(1:100000, 1) #62051
set.seed(62051)
# Read in shapefiles from teams data folder
sara_bound1 <- st_read("./data/SARA_boundary_4269.shp")
sara_veg1 <- st_read("./data/SARA_Veg.shp")
# Check that the projections match; fix the one that doesn't match
st_crs(sara_bound1) == st_crs(sara_veg1) # FALSE- projections don't match.
# sara_bound1 needs to be reprojected to UTM Zone 18N NAD83.
sara_bound <- st_transform(sara_bound1, crs = 26918)
st_crs(sara_bound) == st_crs(sara_veg1) # TRUE
# Quick plot
plot(sara_bound[1])
plot(sara_veg1[1]) # bigger extent than boundary
# Intersect boundary and veg to be same extend
sara_veg <- st_intersection(sara_veg1, sara_bound)
plot(sara_veg[1])
# View attribute table of layers
head(sara_bound) # 1 feature with 95 fields
str(sara_veg)
head(sara_veg)
names(sara_veg)
table(sara_veg$ANDERSONII)
# Simplify vegetation types for easier plotting
dev <- c("1. Urban or Built-up Land", "11. Residential",
"12. Commercial and Services", "13. Industrial",
"14. Transportation, Communications, and Utilities",
"17. Other Urban or Built-up Land")
crop <- c("21. Cropland and Pasture",
"22. Orchards, Groves, Vineyards, and Nurseries",
"31. Herbaceous Rangeland")
shrubland <- c("32. Shrub and Brush Rangeland")
forest_up <- c("41. Deciduous Forest Land", "42. Evergreen Forest Land",
"43. Mixed Forest Land")
forest_wet <- c("61. Forested Wetland")
open_wet <- c("62. Nonforested wetland", "62. Nonforested Wetland")
water <- c("5. Water", "51. Streams and Canals", "53. Reservoirs")
unk <- "Unclassified"
# Create 2 fields in the veg attribute table: simp_veg, and fills
sara_veg <- sara_veg %>%
mutate(simp_veg = case_when(ANDERSONII %in% dev ~ "Developed",
ANDERSONII %in% crop ~ "Open field",
ANDERSONII %in% shrubland ~ "Shrublands",
ANDERSONII %in% forest_up ~ "Forest",
ANDERSONII %in% forest_wet ~ "Forested wetland",
ANDERSONII %in% open_wet ~ "Open wetland",
ANDERSONII %in% water ~ "Open water",
ANDERSONII %in% unk ~ "Unclassified",
TRUE ~ "Unknown"),
fill_col = case_when(simp_veg == "Developed" ~ "#D8D8D8",
simp_veg == "Open field" ~ "#f5f0b0",
simp_veg == "Shrublands" ~ "#F29839",
simp_veg == "Powerline ROW" ~ "#F9421D",
simp_veg == "Forest" ~ "#55785e",
simp_veg == "Forested wetland" ~ "#9577a6",
simp_veg == "Open wetland" ~ "#c497d4",
simp_veg == "Open water" ~ "#AFD0F2",
simp_veg == "Unclassified" ~ "#ffffff"))Sidebar on simple features and shapefiles
Before moving to the next step, I wanted to show how easy it is to create simple features from dataframes that contain X/Y coordinates. We’ll read in a fake dataset in the GitHub repo for this training, and call it wetdat. The dataset contains fake species composition data for wetland monitoring sites in ACAD and includes X and Y coordinates. We’ll use dplyr functions to calculate the number of invasive and protected species in each site by year combination. Then we’ll make it a simple feature, and save it as a shapefile. Note that there are multiple formats you can save simple features as. I show the shapefile version, in case you find yourself going between R and ArcGIS.
library(dplyr)
# Import data
wetdat <- read.csv(
"https://raw.githubusercontent.com/KateMMiller/IMD_R_Training_Advanced/main/data/ACAD_wetland_data_clean.csv")
# Summarize so that each site x year combination has 1 row
wetsum <- wetdat %>% group_by(Site_Name, Year, X_Coord, Y_Coord) %>%
summarize(num_inv = sum(Invasive), num_prot = sum(Protected),
.groups = 'drop')
# Check first 6 rows of output
head(wetsum)
# Convert dataframe to sf
wetsum_sf <- st_as_sf(wetsum, coords = c("X_Coord", "Y_Coord"), crs = 26919)
# ACAD is UTM Zone 19N NAD83, hence the difference between SARA, which is 18N.
# Write sf to disk as shapefile
st_write(wetsum_sf, "ACAD_wetland_data.shp")Step 2. Generate GRTS Sample
The spsurvey package has been updated to point to sf instead of rgdal. It’s a code-breaking change if you have old R scripts that generated GRTS samples. However, the process is even easier now, as you can see below. Here we’re just going to use the simplest of GRTS designs. The output from grts() has multiple slots. The one we want is sites_base, and you can see how we get that slot in the code below.
Note: While I followed the approach documented in the spsurvey vignette to generate reproducable GRTS samples, it does not appear to be working as I’m testing it. Despite running set.seed(62051) after loading spsurvey and then running the code chunk below, each sample appears to be different.
sara_grts <- grts(sara_bound, n_base = 100) # generate 100 points within SARA boundary
sara_grts$sites_base$priority <- as.numeric(row.names(sara_grts$sites_base)) # add priority number (same as row.name)Step 3. Spatial join and filter
First step here is to spatially join the sara_grts layer to the sara_veg layer. Here we only cared about the sara_veg’s simp_veg field, so I used dplyr’s select() function. After joining, you should see a simp_veg field added to the end of the grts layer (it’s actually to the left of geometry, which is the last). In the next step, we then filter the points in the newly created grts_veg layer to only include points that fall in forest habitat.
# Spatial join
grts_veg <- st_join(sara_grts$sites_base, sara_veg %>% select(simp_veg))
names(grts_veg)
# Create list of forest habitats
sort(unique(sara_veg$simp_veg))
forest_veg <- c("Forest", "Forested Wetland")
# Filter out non-forest points
grts_forest <- grts_veg %>% filter(simp_veg %in% forest_veg)
nrow(grts_veg) # 100 points
nrow(grts_forest) # fewer pointsStep 4. Create map
Now it’s time to plot. The code below may seem a bit much, but we’ll break it down piece by piece. The great thing about building maps this way is that you’re essentially building a template in code that you can steal from and adapt for future maps. You can even turn these into a function that’s even easier to use in future maps. That’s beyond what we can cover here, but it’s another great benefit of making maps in R instead of ArcGIS.
The code below is broken into the various pieces that make up the map. The way tmap works is that first, you have to add the layer via tm_shape(), and then you specify how you want that layer to look, like tm_fill(), or tm_border(). Each piece has its own legend as well. This is similar to how you start each ggplot2 graph defining the data with ggplot(data, aes(x = xvar, y = yvar)), and then start adding geom_point(), or whatever else to define how the data are plotted. The difference with tmap is that every layer you want in the map has to be coded this way. Finally tm_layout() is similar to theme() in ggplot2, and is where you can customize the map layout.
-
The first line that defines
for_legendmakes a list of the habitat types in simp_veg and their associated colors. That saves time having to type them all over again, and potentially spelling them wrong. -
The first
tm_shape()in the map sets the projection and the bounding box. If you don’t set the bounding box, the map will show the largest extent of your layers. So if you have a roads layer at the state-level, your map will be zoomed at that extent, instead of the park boundary. -
Note the difference between
tm_fill(), which fills in colors, andtm_borders(), which only plots outlines. If you want both borders and fills, usetm_polygons()instead. -
The z value in the
tm_add_legend()allows you to set the order that each legend appears in the map. Otherwise, they’ll appear in the order you specify the layers. -
The code for
tm_symbols()allows you to change the symbol format in a similar way toggplot2. We then addedtm_text()to label each point using the numbers in the priority column of the grts_forest attribute table. Thexmodandymodallow you to offset labels from the points either horizontally and vertically. In this case, negativexmodmoves the label to the left, and a negativeymodmoves the label below the point. The default location for labels is directly on top of the point. -
The
tm_layout()code is where you can change the default settings of the figure, including font size, placement of the legend, and margins. The title name and position are also set in thetm_layout(). -
The
tmap_save()allows you to write the map to disk and to specify the height and width of the map. I prefer to write maps to disk to see what the size looks like before changing any of the layout and font size settings, because the figure on your disk will look different (and often better) than the plot view in your R Studio session.
# Creating list of simp_veg types and their fill colors for easier legend
for_legend <- unique(data.frame(simp_veg = sara_veg$simp_veg, fill_col = sara_veg$fill_col))
sara_map <-
# Vegetation map
tm_shape(sara_veg, projection = 26918, bbox = sara_bound) +
tm_fill("fill_col") +
tm_add_legend(type = 'fill', labels = for_legend$simp_veg,
col = for_legend$fill_col, z = 3) +
# Park boundary
tm_shape(sara_bound) +
tm_borders('black', lwd = 2) +
tm_add_legend(type = 'line', labels = "Park Boundary", col = "black",
z = 2)+
# GRTS points
tm_shape(grts_forest) +
tm_symbols(shapes = 21, col = "#EAFF16", border.lwd = 0.5, size = 0.3) +
tm_text("priority", size = 0.9, xmod = -0.4, ymod = -0.4) +
tm_add_legend(type = 'symbol', labels = "GRTS points", shape = 21,
col = "#EAFF16", border.lwd = 1, size = 0.5,
z = 1) +
# Other map features
tm_compass(size = 2, type = 'arrow',
text.size = 1, position = c('left', 'bottom')) +
tm_scale_bar(text.size = 1.25, position = c("center", "bottom")) +
tm_layout(inner.margins = c(0.2, 0.02, 0.02, 0.02), # make room for legend
outer.margins = 0,
legend.text.size = 1.25,
legend.just = 'right',
legend.position = c("right", "bottom"),
title = "Saratoga NHP GRTS points",
title.position = c("center", "top"))
tmap_save(sara_map, "SARA_GRTS.png", height = 10.5, width = 8,
units = 'in', dpi = 600, outer.margins = 0)
Here’s what sara_map looks like in your plot viewer:
sara_map
Here’s what the map looks like after it’s written to disk and the dimensions are set using tmap_save(): 
The last cool thing to show with tmap, is that you can include interactive HTML widgets similar to what leaflet can do (coming next). With the interactive mode, you can change baselayers, turn your layers off/on, and zoom to different extent. The legend in the interactive mode is a bit limited to only show fills, but it’s still pretty cool.
tmap_mode("view") # turn interactive mode on
sara_map